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ABSTRACT 


^ 

First  an  investigation  of  modeling  stochastic  processes 

by  difference  equations  (Markov  process)  was  undertaken. 

The  starting  point  of  the  modeling  procedure  is  the  knowledge 

of  the  spectrum  of  the  process.  Two  methods  are  discussed. 

One  is  based  on  optimal  estimation  theory  and  leads  in  most 

cases  to  a high-order  (perhaps  infinite)  Markov  process.  The 

second  method,  based  on  linear  system  theory,  leads  to  a 

first  order  Markov  process  (in  matrix  representation) . Both 

methods  have  been  extended  to  two-dimensional  processes. 

Secondly,  recursive  estimation  (filtering)  of  two-dimensional 

random  fields  was  addressed.  It  was  shown  that  a t.>ro- 

dimensional  recursive  filter  cannot  be  optimal.  Therefore, 

only  a sub-optimal  solution  is  available.  This  solution 

minimizes  the  mean  square  error  for  a specific  structure  of 

a filter.  Finally,  applications  of  modeling  and  recursive 

filtering  are  discussed.  An  image  that  includes  a target, 

correlated  noise  and  random  noise  was  processed.  Some 

methods  of  target  enhancement  (also  called  "restorati^on" 

are  discussed. 
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NOTATION 


x(k,A) 

u(k,Jl) 

v(k,i.) 

X (k, H) 
y{k,l) 

T 

Og,a 

Pl»p2 

®1'®2 

TH 

Pg(k,Jl) 

P(k,i) 

E(k,a) 

a 

n 

R 

o 

u 

Q 

P.  .(k,^) 
/ J 

(.)^ 

X 

G 


correlated  field 

modeling  error  (random  wlak  of  the 
Markov  process) 

white  noise  field. 

estimation  of  x(k,6) 

^ ^x(k,l)  - x(k,£)]  = estimation  error 

A (ftk+l,Z)  ^(k,Z)  £(k,Jl+l)) 

target  intensity 

variance  of  correlated  process 

correlation  coefficients  that  determine 
the  bandwidth  of  the  process 

correlation  coefficients  that  determine 
the  center  frequency  of  the  process 

value  of  threshold  in  the  decision  process 

- variance  of  error 

2 

» (k,Z)  = meam  of  square  erros 

A E{^(k,£)-|^{k,J?.)  } 

variance  of  random  noise 

a 2 
n 

- variance  of  modeling  error  (random  walk) 

A a 2 
u 

A E{(?(k,)l)*  <?(k+i,  l+j)  } 

- means  the  transpose  of  the  vector  (•) 
mean  of  correlated  field 

- steady  state  gain 


_ . . . a. 


7 


I.  INTRODUCTION 


1.  The  goal  in  this  thesis  was  to  solve  the  problem 
described  in  the  following  paragraphs.  Consider  an  image 
sensing  device  that  has  an  array  of  NxN  sensing  elements. 

The  images  that  are  sensed  by  the  device  include  some 
targets  that  suffer  from  degradation  due  to  background 
(for  exaunple,  clouds).  The  background  noise  eliminates  the 
possibility  of  direct  detection  of  the  target.  The  background 
is  assumed  to  be  a correlated  noise  source. 

Further  interference  to  the  output  of  the  imager  might 
be  the  internal  noise  of  the  device.  It  is  assumed  that 
the  noise  of  each  sensing  element  is  uncorrelated  to  the 
others  (white  noise).  Therefore,  the  output  of  the  imager, 
y(k,£.)  includes  three  types  of  processes: 

1.  the  target  T(k,£) 

2.  The  correlated  background  x(k,)l) 

3.  The  white  noise  v(k,i) 

y(k,Z)  = T(k,)l)  + x(k,«,)  + v(k,il)  . 

A typical  image  is  seen  in  Figure  2.  The  problem  in  this 
thesis  is  to  detect  the  intensity  and  location  of  the  target 
by  using  recursive  techniques  that  are  applicable  for 
real-time  hardware. 


8 


2.  The  solution  to  this  problem  is  carried  out  in  three 
steps : 


a)  Mathematical  modeling  of  the  background. 

b)  Estimation  of  the  background. 

c)  Detection  of  the  target. 

The  main  idea  was  to  eliminate  the  background  by  subtraction 
of  the  estimated  background  from  the  original  imager's  output. 
Then  the  residual  image  includes  only  the  target  and  the 
white  noise.  The  detection  is,  therefore,  easier.  Still, 
for  detection  with  small  errors  (false  alarms  and  misses) 
the  target  must  have  an  intensity  higher  than  the  white 
noise.  The  detection  procedure  is  shown  in  Fig.  1. 


Fig.  1:  Target  Detection  from  a Noisy  Image 

In  each  of  the  three  steps,  some  original  ideas  have  been 
developed. 
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3 . In  the  "mathematical  modeling  of  the  background"  the 
starting  point  is  that  there  is  a correlation  between  points 
of  the  background.  The  two  dimensional  autocorrelation 
function  was  assumed  to  be  known.  Only  stationary  random 
fields  were  treated. 

The  main  idea  of  the  mathematical  modeling  is  to  repre- 
sent the  image  by  difference  equations  that  are  forced  by 
white  noise.  The  conventional  way  of  representing  a process 
whose  autocorrelation  function  is  given  is  by  using  a 
"Markov  Process."  The  value  of  a point,  x(k,2.),  is 
represented  by  its  neighbors: 


x(k,Jl)  = Za  x(p,q) 

P/  4 

{(p,q)}  is  a group  of  "neighbors"  near  the  point  (k,il).  And: 

x(k,£)  = x(k,i!,)  + u(k,Jl) 

=.  ZOp  gX  (p,q)  + u(k, «,)  . 

u(k,H)  is  called  "modeling  error." 

The  weighting  coefficients  a are  chosen  to  make  the 
variance  of  the  modeling  error  minimum.  This  is  done  by 
using  the  well  known  "orthogonality  principle"  in  optimal 
estimation  theory  [see  ch.  II  section  H] . This  technique 
suffers  from  some  difficulties  as  follows: 
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have  no  simple 


The  weighting  coefficients  a 
expressions. 

The  number  of  neighbors  one  has  to  use  is  theoretically 
infinite  for  some  types  of  processes. 

- The  coefficients  are  found  by  solving  many  algebraic 
equations,  especially  in  the  two-dimensional  case. 

The  number  of  equations  is  equal  to  the  number  of 
coefficients . 

It  is  difficult  to  describe  the  two-dimensional 
difference  equation  in  a "state-vector  structure". 

The  optimal-estimation  approach  to  linear  modeling  is 

summarized  in  Ch.  II  section  H. 

The  method  of  modeling  that  was  derived  in  this  thesis 

is  to  find  a recursive,  linear,  invariant  filter 

such  as: 

I - when  forcing  the  input  with  white  noise,  the  output 

of  the  filter  will  be  a random  process  that  has  the 
same  autocorrelation  function  as  the  given  field. 

So,  instead  of  dealing  with  the  spectrum  of  the  process, 

one  has  to  deal  with  white  noise  and  a linear  filter.  This 

is  well  known. 

Rosser  [Ref.  12],  in  1975,  suggested  a state-space- 
structure  for  two  dimensional  fields  and  also  derived  some 
of  the  properties  for  this  structure  [see  Ch.  II  Section  G] . 

All  of  the  examples  that  were  chosen  in  this  thesis  led  to 
a filter,  ¥l(z^,Z2)  , whose  state-space-structure  fitted 
easily  into  Rosser's  structure. 

Moreover,  this  method  of  modeling  doesn't  suffer  from 
any  of  the  disadvantages  of  the  previous  method.  It  should 
be  emphasized  that  this  method  is  limited  only  to  separable, 
two-dimensional  autocorrelation  functions  and  to  causal  models. 
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This  method  is  called  "Filter  Response  Method"  and 
is  summarized  in  Ch.  II  Section  F. 

4.  The  estimation  of  the  background  is  actually  an  exten- 
sion of  the  one  dimensional  Kalman  Filter.  The  problem  is 
roughly  defined:  "Given  a noisy  measurement 

y(k,£)  = x(k,Jl)  + v(k,)l) 

where  x(k,)l)  is  a correlated  field  and  v(k,£)  is  white 
noise.  The  filter  should  estimate  x(k,2.),  denoted  x(k,£). 

It  was  shown  that  a recursive,  two  dimensional  filter  cannot 
be  optimal  in  the  sense  that  the  error  cannot  be  made 
orthogonal  to  the  measurements,  as  some  researchers  have 
tried  to  do  [Ref.  10].  There  is  a distinct  difference 
between  one  dimensional  processing  and  two  dimensional 
processing.  Therefore  the  method  in  this  thesis  is  to 
define  a reasonable  structure  for  the  recursive  filter,  and 
then  to  calculate  the  parameters  of  the  filter  to  minimize 
the  mean  of  the  square  error.  Recursive  equations  for 
calculating  the  filter  parameters  (gain)  and  variance  of 
error  were  developed. 

The  results  were  checked  in  two  ways : 

1)  By  comparing  with  the  optimal,  non  recursive 
estimator. 

2)  By  simulation.  A correlated  image  was  added 

to  a field  of  white  noise.  Then  the  correlated 
part  of  the  combined  image  was  estimated  by 
using  the  recursive  filter.  The  variance  of  error 
was  compared  with  the  theoretical  variance  of  error. 

The  results  for  all  cases  that  were  checked  show  good  coincidence. 
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< 5.  In  the  study  of  target  detection,  the  problem  was  to 

detect  lines.  Most  of  the  algorithms  suffer  from  one 
disadvantage.  When  detecting  a line,  an  a-priori  assumption 
must  be  made  about  the  direction  of  the  line.  Therefore, 
in  order  to  detect  lines  in  several  directions,  the  image 
has  to  be  scanned  several  times,  or  during  one  scanning  to 
do  several  calculations  for  each  point.  No  doubt  that  for 
real  time  applications  all  of  those  algorithms  have  a great 
disadvantage.  The  algorithm  developed  in  this  thesis 
detects  lines,  regardless  of  their  direction.  The  key 
point  in  this  algorithm  is  feedback  from  the  detection 
(which  is  a decision  process)  to  the  filter.  Most  algorithms 
[ of  filtering  and  detection  do  it  in  two  separated  steps. 

Fig.  3 shows  the  target  that  was  detected  from  the  original 
image  in  Fig.  2. 

6 . Ou -line  Of  Chapters 

Chapter  II  discusses  the  problem  of  two  dimensional 
random  fields.  In  order  to  make  this  chapter  "stand  alone" 
for  reading,  some  background  material  was  included  in  Sections 
I B and  C.  This  background  material  includes  information  on 

two-dimensional  operations  and  random  processes.  Section  E 
describes  models  for  a one-dimensional  random  process. 

Section  F extends  the  methods  of  modeling  for  two  dimensional 
processes. 

Chapter  III  is  a review  of  estimation  theory.  It  explains 
the  "position"  of  recursive,  linear  estimation,  with  quadratic 
form,  in  estimation  theory. 
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Chapter  IV  solves  the  problem  of  estimating  a two- 
dimensional  process,  when,  originally,  this  process  is 
combined  with  white  noise. 

Chapter  V shows  the  application  of  Ch.  II  and  IV 
(modeling  and  estimation)  to  the  specific  problem  of 
target  detection. 
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Fig.  2:  A Typical  Image,  Includes  a Target, 

Correlated  Noise  (Background)  and 
Random  Noise 


II.  LINEAR  MODELING  OF  TWO-DIMENSIONAL 
RANDOM  FIELDS 


A.  INTRODUCTION 

1.  This  chapter  contains  the  basic  principles  of  modeling 
a random  process  by  linear  equations.  The  basic  assumption 
will  be  that  the  statistical  relationship  between  each  two 
measurements  is  known.  This  relationship  is  the  so-called 
autocorrelation- function.  The  techniques  that  will  be  used 
are  mainly  extensions  of  the  principles  that  are  used  in  the 
one-dimensional  case.  Terms  like  "Markov  Model",  "state 
variables",  "Factorization"  will  be  used. 

2.  In  the  one -dimensional  case  the  state-equations  have  the 
foirm: 


x(k+l)  = A*x(k)  + B*u(k) 

The  words  "one-dimensional"  refer  to  the  variable  k (but  the 
state  vector  x(k)  can  be  a multidimensional  vector) . It  will 
be  shown  how  to  use  state-equations  for  a two  dimensional 
random  field,  in  which  the  variables  are  k and  il  . 

Note ; When  a random  variable  depends  only  on  one  dimension, 
k,  it  is  called  a "stochastic  process".  When  the 
dimensionality  is  two  or  more,  as  k,£,  it  is  called 
a "random  field". 
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3.  What  is  the  philosophy  behind  describing  a random  process 
by  state  equations?  The  answer;  it  is  difficult  to  deal  with 
random  fields  in  which  there  is  correlation  between  adjacent 
points.  Therefore,  the  main  idea  is  to  show  how  a correlated 
field  can  be  described  as  an  output  of  a_  linear  system  which 
is  driven  by  white  noise.  So,  instead  of  handling  a 
correlated  field,  one  has  to  hcindle  white  noise  and  a linear 
system.  The  last  problem  is  well  known  from  system  theory. 
Figure  4 illustrates  the  concept  of  modeling  by  the  "filter 
response  method" . 


INPUT  IS 
WHITE  NOISE 


LINEAR 

filter 


OUTPUT  IS  THE  GIVEN 
CORRELATED  FIELD 


Fig.  4:  Concept  of  "Modeling"  a Random  Field  as  an 

Output  of  a Linear-Filter 


4.  Although  this  thesis  mainly  delas  with  discrete  fields, 
continuous  models  will  be  described  also,  in  some  cases. 

5.  In  order  to  make  this  chapter  "stand  alone"  for  reading, 
some  background  material  was  included  in  sections  B,  C. 

6.  Section  D describes  the  "Filter  Response  Method"  of 
modeling  a random  process. 
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7.  Section  E describes  the  one  dimensional  processes,  and 
Section  F is  its  extension  to  two  dimensions.  The  Markov- 


Process  is  the  key  to  the  modeling  procedure.  It  will  be 
shown  that  from  the  same  starting  point  (the  given  auto- 
correlation function)  two  models  can  be  developed. 

8.  The  modeling  by  the  "Filter  Response  Method"  leads  to 
models  that  have  the  structure 


• 

M(k+l,i) 

A i 

M(k,£) 

- 

N(k,)l+1) 

• 

N(k,^) 

• 

j 

• u(k, i) 


Roesser  [12]  described  the  properties  of  that  structure  and 
section  G is  a summary  of  those  properties. 

9.  Section  H summarizes  another  method  of  linear  modeling, 
using  optimal  estimation  theory  (the  orthogonality  principle) . 
This  method  leads  to  different  results  than  the  filter  response 
method  with  the  exception  of  one  case. 

10.  So,  two  methods  of  modeling  will  be  shown  (sections  F,  H) . 


Section  I is  a summary  of  this  chapter.  It  also  compares 
the  two  modeling  methods. 


B . BACKGROUND  MATERIAL : 

TWO  DIMENSIONAL  OPERATIONS  AND  FILTERS 

This  section  is  a siommary  of  two  dimensional  operations 
as  Fourier  transformation,  2-  transformation,  two  dimensional 
filters,  etc.  Some  of  these  operations  are  used  in  the  next 
sections . 

1.  Fourier  Transformation 

Definition;  given  a function  h(x,y),  the  Fourier 
transformation  H{a)  ,o)  ) is 

X y 


H (u)  , oj  ) = 

X y 


/ / h(x,y)  • e 


• e ^ - dxdy 


= T'  I h(x,y)| 


and  the  inverse  transform 


jo)  y 

h(x,y)  = / /H(oi^,(i)^)*e  *e  ^ ^^y 


The  two  dimensional  Fourier  Transformation  is  mostly  used  in 
spatial  operations,  and  ,01^  are  called  spatial  frequencies , 
The  properties  of  the  Fourier  transformation  are; 


Linearity  ^^a*f^(x,y)  +b*f2(x,y)| 

» t b*F2 
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Scaling 


f (aX^ 


2-4 


Shift  Operation 


f (x-a 


y-b)( 


-j  (ui  a+u)  bj 

e ^ ^ ‘ F (u)  , u)  ) 

X y 


2-5 


Convolution 


i 

I 


f (x,y) *h (x. 


2-6 


Parsaval  Theorem 


00  00 


00  00 


/ / f (x,y) • g* (x,y)  dx  dy 

— oo  -»co 


/ / F (cj  ,U1  ) - G*  (u  ,0)  )’du}  du) 

X y * y ^ y 

*00  *oo 

2-7 


Autocorrelation 


/ / f (C,n) -f  (5-x,n-y), 


*00  *00 


2 

lF(u)^,ujy)|  2-8 
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Gradient 


X X y 


r 


Inversion 


■i 


L 


.y)| 


= r 


r(f<x, 


y) 


= f(x,y) 


Rotation 


f (x,y)V  [ = f(-x,-y) 


2-9 


2-10 


2-11 


2-12 


Equation  2-9  will  be  useful  in  the  next  sections  and  the 
proof  for  this  equation  is  given  as  follows: 


r(^ 


/ / ^_f(x,y)  . 

iX 


X -DU  y 

e • e ^ dxdy 


.oo  «.oo 


I 


00  r~  00 

/I  f . 


ox 


e * dx 


-3“„y 


• e 


dy 


— 00  [—00 
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Using  integration  by  parts: 


3F(x,y)/  / 

— '=  / f(x,y).e  "" 


+ jw  • / f (x,y)  • e 

X 


-ICD  X 
X 


-JO)  y 

■e  y • dy 


If  f(x,y)  = 0atx=±'»  ^ y = 


ju)^-  / / f (x,y)  *e 


-jw  X y 


= j (Jj  • F (d)  , OJ  ) 

X X y 


Note:  In  the  one  dimensional  case  most  equations  appear  with 

"t"  as  the  variable.  In  order  to  make  the  equations 
in  the  two  dimensional  case  similar  to  the  one 


dimensional  case,  the  notations  x,  y,  oj  , oj  will 

X y 

be  replaced  by  t^^,  t2t  0)2  • 


i 
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2.  Two  Dimensional  Z Transformation 


Definition:  Given  a discrete  field  X(nj^,n2),  the 

two  dimensional  z transformation  is  defined: 


CO  CO 


^{x(n^,n^)}  = X(z^,z^)  = ^ ^ x (n^,n2> • 

2- 


n^=-oo  n^=-<^ 


An  important  property  of  the  transformation  which  will  be 
used  later  is 


if;  -2”[x(n^,n2)] 


X(z^,Z2> 


2- 


then : 


fx(nj^-l,n2)| 


Zj^”^X(Zj^,Z2) 


2- 


2"  |x(nj^,n2-l)j 


Z2"‘^X(z^,Z2) 


2- 


Proof : {x(n^,n2)  } 


-n  -n2 

I I x(nj^-l,n2)  Z^  -Z2 


«.00  ^00 


z I X(m  ,n2>  • z^~  • Z2  ^ 


mm  CO  — QO 


OO  00 


-1  ""2 

Zj^  EE  x(m  ,n2)  • Z2 


-•00  *00 
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i ^ ' -B 


-n2 

2 

13 

14 

15 

16 


• f f x{n^,n^) 

«_oo  «oo 


= z^^.%  {x{n^,n^)] 


3.  Pasponse  of  a Linear  System 

If  h(tj^,t2)  is  the  "point  spread  function"  of  a two 
dimensional  system,  then  the  response  to  a given  input  x(t^,t2> 
will  be  the  two  dimensional  convolution: 


y(t^,t2)  = / / h(tj^-T^,t2~T2)  *x(t^,t2)  *<iT2dTj^  2-17 


_00  .00 


and  the  transfer  function  is  defined: 


H (ojj^,aj2) 


Y ((jjj^,o)2) 
X ((jj^,o)2) 


2-18a 


In  the  discrete  case,  the  convolution  has  the  form: 


00  00 


y(k,£)  = E Z h (n  ,m)  • X (k-n , ii-m) 


2-19 


m=-oo  n=-« 


where  h(n,m)  is  the  discrete  "point  spread  function"  and 
the  discrete  transfer  function: 

Y(z. ,z_) 

»(^1'^2>  • xri:75J  2-18b 


•1'  2' 
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4 . Separable  Functions 

Definition;  f(tj^,t2)  j-s  called  "separable"  if 
are  two  one-dimensional  functions  f^(t),  ^2^^^  such  as 

f(ti,t2)  = f (t^)  • f 2 (t2) 

Theorem 


If:  f(tj^,t2)  = f (t^)  • f 2 (t2) 

Then:  F(uJj^,(jJ2)  = (cj^)  • F2  (^2 ) 

00 

— j (1)  t 

where:  F^(o)j^)  = / fj^(tj^)*e  ^ ^‘^^1 


-e  ^ ^.dt2 


In  the  discrete  case: 


If;  x(n^,n2)  = X (n^^)  • X2  (n2) 


Then;  2^{x(n^,n2))  = -2’{x  (n^^)  } *^{x  (n2>  } 


there 


2-20 


2-21 
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C.  BACKGROUND  MATERIAL: 

TWO  DIMENSIONAL  LINEAR  SYSTEMS  WITH  RANDOM  INPUTS 
In  this  section  expressions  for  the  response  of  two 
dimensional  filters  to  random  inputs  are  developed.  It  will 
be  shown  that  the  expressions  are  very  similar  to  the  one 
dimensional  case.  First,  a brief  review  of  analysis  of 
random  signals  (in  two  dimensions)  is  introduced. 

Given : A "brightness  function"  of  a two  dimensional  field, 

X (tj^,  t2>  [see  fig.  5]  . 


Fig.  5:  The  Two  Dimensional  Field 

1.  The  probability  that  x(tj^,t2)  lies  in  a finite  range 
a £ x(tj^,t2)  £ b at  a point  in  the  field,  (tj^,t2),  is  given 
by  the  integral 

b 

Pr(a  £ x(tj^,t2)  £ b)  » / f (x,  tj^,t2)  • dx  2-22 

a 
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► 

I 

I 


fl(x,t^,t2)  is  the  "probability  density  function." 

2.  The  first  moment  is  the  expected  value  (or  mean)  and 
is  defined  by  the  integral: 

CO 

xTE^Tt^  = E{x(t^,t2).}  = / X*  f (x,  tj^,  t2)  *dx 

•—00 

2-23 


3.  The  second  order  moment  (variance)  of  a point  in  the 
field  is: 

E{x  (tj|^,t2)}  = X (t^,t2)  = / x^*  f (x,t^,  t2)  *dx 

••00 

[ 2-24 

4.  The  autocorrelation  function,  between  points  (t^^,t2)  and 
(ti',t2')  is  defined  as: 

^xx^^l'^2'^l' '^2*  ^ x(tj^,t2)  • x(tj^\t2’  ) 

00  00 

= J j"  /f2f^l^tj^,t2),Xj^’  (t^' ,t2*  ))  dxdx’  2-25 


J 
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1 


i 

r 


Fig.  6;  The  Parauneters  That  Take  Part  in 
the  Autocorrelation  Equation 


£2^’)  is  the  "joint  probability  density  function". 

Fig.  6 shows  the  parameters  that  take  part  in  the 
autocorrelation  equation. 

5.  Stationary  Fields. 

The  assumption  that  the  field  is  stationary  means  that 
the  statistics  of  a point  in  the  field  is  not  dependent 
on  the  location  of  the  point. 

In  this  case  the  mean  and  variance  have  the  form 

E{x{tj^,t2)}  * X ^ Etx^  (tj^,t2)  } = 2-26 

and  the  autocorrelation  function; 


1»T2- 


00  00 

/ / X • x'  ’f2(x  ,x* , X2)dx*dx'  2-27 

wOO  *00 
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Such  a field  is  called  (in  the  two  dimensional  case)  a 
"Homogeneous  Field." 

6.  Ergodicity 

A further  simplifying  assumption  which  is  usually  done 
during  the  analysis  of  stationary  signals  is  the  ergodic 
property.  This  hypothesis  states  that  under  certain  con- 
ditions, present  in  many  cases,  the  statistical  averaging 
of  X at  a given  point  is  equal  to  the  spatial  averaging 
of  all  points.  That  means: 


= - - 


<V 


2-28 


where  by  definition: 


^ f f x(t  ,t2) -dtj^*dt2  2-29 
2 — » 


<^x)  i 


is  the  spatial  averaging. 

The  autocorrelation  function  will  be: 


T2) 


E{x(tj^,t2)  X ( tj^+Tj^,t2+T2)} 


2-30 


lim 


T7t  Tt~ 


f f x(t^,^2>  X (tj^+Tj^,t2+'^2^'^^l‘^^2 


■ 2-^* 
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^ . Power  Spectrum  Density  Function 

In  the  study  of  random  signals  the  concept  of  power 
spectral  density  function  takes  place.  For  the  purpose  of 
this  review  the  spectral  density  of  two  dimensional  fields 
is  defined  as  the  two  dimensional  Fourier  transof rmation 
of  the  autocorrelation  function.  It  is  called  P (cuwW-)  : 

XX  X X 


'■’“<’'xx<"l'''2>> 


= / / Rxx^^1''^2^ ® -dT^.dT^  2-31 

« 00  00 

and  the  inverse  transform  gives  the  result: 

® « j(jj 

R (T  ,T_)  = f f P (u),,oj_)»e  *e  •du)T*dcu- 

XX  12  „_xx  l 2 1 2 

wOO  »00 

2-32 

From  the  last  expression,  P ,w_)  gives  the  density  of 

mean  square  value  of  the  variable  over  the  spectrxim  of  the 
real  frequencies. 

8.  Response  Of  A Linear  System  To  Random  Inputs 

Given:  A linear,  two-dimensional  linear  filter 

with  a transfer  function  H ( ju)2)  • The  input  to  the 
filter  is  a stationary  process  x(tj^,t2)  with  an  autocorrelation 
function  R 

XX  X V 
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Fig.  7:  Response  Of  A Filter  To  Random  Input 

Question:  Find  the  autocorrelation  of  the  output  signal. 

The  answer  will  be  presented  in  n6xt  theorems. 


Theorem ; The  mean  of  the  output  is  given  by 


QO  00 

= / Z (tj^-a^,t2-a2)  ^-h(aj^,a2) -da^*da2 

and  because  the  process  is  stationary: 


E{x(tj^-aj^,t2-a2)  } » E{x(tj^,t2)}  = x 
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y = / f K ' h(a,,a_)  • da, "da- 

.00  .CO  X A X ^ 


X f f h(a, ,a») *da,  *da- 

-CO  _o>  1 ^ 1 Z 


X f f h(a,/a2)*e  • e • da,  -da- 

— 00  —CO  ± ^ X ^ 


H(  joj2) 


0)^=0 

u)2=0 


Q.E.D. 


Theorem;  The  cross-spectra  of  x,y  is  given  by 


P^y(0ii,<U2)  = * H(ja.^,joj2) 


“ E{yCt^+Tj^,t2+T2)  *x(tj^,t2)  } 


E(y(t^,t2)  * x(t^“Tj^,  t2“T2)  } 
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• _^-r> 


00  00 


f 


f E{x(t^-T^,t2-T2) •X(t^-a^,t2~a2) | *h(a 


R^y(Tl,T2)  = 


XX 


(t^, T2) 


h(Ti,T2) 


The  last  expression  is  a convolution.  Therefore, 
the  Fourier  transform  leads  to; 


^xy^‘^l'‘"2^  = ^xx^‘“l'‘^2^  * H(ju)^,ja)2) 

Q • £ • 0 • 


a2)da^da2> 


,a2)  duj^da^} 


,a2)da^da2> 


2-35 


taking 


Theorem;  The  spectral  density  of  the  output  is  given  by; 


Proof ; 


= E{y(tj^,t2)-y(t^+Tj^,t2+T2)  } ' 

00  OO 

=■■  E{y(tj^+T^,t2+T2)  ' / / x(t^-aj^,t2-a2)  •h(aj^,a2)  *da^da2} 

“•OO  “OO 

00  OO 

= / / R (a^+f^,a2+'^2^ 

W30  «^0  * 

define:  = -62 

da^  = -d6j^  da2  = -d02 


00  00 


/ d9j_de. 


Ryy(Ti,T2)  = Ryx<"l'"2^  ^ h(-T;^,-T2) 


2-37 


Taking  the  Fourier  transform  of  2-37; 


Pyy<‘^l'‘^2^  * ^yx^‘^l'‘^2^  -H  (-ju)j^,-ju)2) 


But; 


Py^(u)i,a)2)  . P^^  • H(ju.^,ju)2) 
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and : 


Pyy(aj^,a)2)  = ' « ( joj^,  j(ij2)  • H (- jco^ , > jco2 ) 

= ‘H  (ju)^,  joj2)  *H*  (jo)^,  jo)2) 

2 

= ^xx^‘^l'‘^2^  ’ H jo)2) 

Q • E • D • 

The  Discrete  Case; 

In  the  discrete  case,  a stationary,  two-dimensional 
field  has  an  autocorrelation  function  defined  as: 

R (n,m)  * E{  X (k,  «,)  X (k+n,  Jl+m)  } 2-38 

XX 

Equivalent  to  the  power  density  function  in  the  continuous 
case,  we  define: 

Since  it  is  known  that  the  z transformation  is  related  to 
the  Fourier  transform  by: 


36 


z 


- j • (!)•  T 
e 


It  is,  therefora,  clear  that  after  passing  a random  process 
through  a two-dimensional  discrete  filter,  the  statistics 
of  the  output  are: 


H(z. fZ,! 


2-40 


2-41 

The  last  two  expressions  can  be  proven  exactly  in  the  same 
way  as  in  the  continuous  case,  by  starting  from  the 
discrete  convolution. 

9.  Isotropic  Fields 

The  homogeneous  (stationary)  field  was  defined  as 
the  case  where  the  statistics  of  the  field  do  not  depend 
on  the  location.  If  the  statistics  are  not  only  independent 
of  the  location,  but  also  independent  of  the  direction , 
the  field  is  called  isotropic. 


i 

_L 
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Example : 


R (t,  ,To)  = e ^ . e 

XX  1 2 


the  field  is  not  isotropic.  The  correlation  is  greater  in 
the  directions  of  the  system  axis,  than  in  other  directions. 
But  if: 


^xx^'^l''^2^  ® 


- "hi 


2~  T 
^^2 


the  field  is  isotropic.  The  correlation  between  two  points 
depends  only  on  the  distance  between  these  two  points  and 
does  not  depend  on  the  direction. 


BACKGROUND  MATERIAL:  THE  CONCEPT  OF  MODELING  BY 

"FILTER  RESPONSE  METHOD"  AND  THE  FACTORIZATION  PROBLEM 

1.  Continuous,  One  Dimensional  Case 

Given:  A covariance  function  of  a process,  R„,^(t,T). 

Find:  A linear  system  (differential  equation  model) 

such  as: 

when  the  input  is  white  noise 

the  out  -lut  has  an  autocorrelation  function 


Fig.  8 defines  the  problem: 
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WHITE  NOISE 


A linear  system 

X = A X + 8U 
IN  FREQUENCY  DOMAIN 
H ( j o)  ) ~ 7 


R^^  (t,  r)  18  given 


Fig.  8:  Definition  of  the  Modeling  Problem 

in  the  Continuous  Case 


Solution : 


First,  the  solution  exists  only  if  x(t)  is 


stationary , say: 


R (t,T)  = R (t: 

XX  ' XX 


In  this  case,  by  using  Eq.  2-36  and  by  assuming  that  the 


transfer  function  of  the  required  filter  is  H(jw),  the  power 


spectral  density  of  the  output  is: 


Pxx^“)  = Pqu  (a))  • H(jtj)  -H*(jw) 


u(t)  is  assumed  to  be  white  noise.  therefore; 


P (oj)  = 1 

uu 


2-42 


= H ( joj)  *H*  ( joj)  = 


given  power 
spectrum 


2-43 


In  this  problem  P..„  (oj)  is  given.  Therefore,  the  solution 

xx 

for  the  required  filter  is  to  find  a function  H(ja))  that 
satisfies  Equation  2-43. 

Such  a solution  exists  for  symmetric  autocorrelation 
functions  (then  P (oj)  is  a real  function  of  cj)  . 

XX 

2 . Discrete,  One  Dimensional  Case. 

In  this  case  R „(n)  is  given. 

XX 

It  is  required:  to  find  a discrete  filter,  H(z), 

so  that  when  the  input  is  white  noise,  the  output  will 

have  the  given  autocorrelation  function  R (n) . 

xy 

Here,  the  solution  is  by  using  the  equation: 

= P^^(z) -H{z)  . 

If  the  input  is  white  noise  then: 


= 1 

Therefore  the  solution  for  the  required  filter  is  to  find 
a function  H(z)  that  satisfies: 

P^^(o))  = H(z)  • H(z'^) 

XX 

where  P_...  (uj)  is  given  in  our  problem. 

XX 
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WHITE  NOISE 


H = ? 
(z) 


( n ) is  given 


Fig.  9:  The  Modeling  Problem  in 

the  Discrete  Case 


2 . The  Two-Dimensional  Case 

In  the  two  dimensional  case  the  extension  of  this 
method  requires  factoring  H 

and  P^^(z^,Z2)  to  H(z^,Z2)‘H(Zj^  ^'^2  problem  is 

that  there  is  no  factorization  technique  in  the  two- 
dimensional  case.  Therefore  this  modeling  method  is 
limited  to  separable  autocorrelation  functions.  Section 
H shows  another  method  of  modeling  (optimal  estimation 
approach)  that  does  not  suffer  from  this  limitation,  but 
has  other  disadvantages. 

E.  BACKGROUND  MATERIAL:  MODELING  OF  ONE  DIMENSIONAL  RANDOM 

PROCESSES  (BY  "FILTER  RESPONSE  METHOD") 

There  is  a special  class  of  stationary  random  processes, 
which  is  very  common.  The  basic  assumption  of  these  random 
processes  is  that  the  correlation  between  two  points  decreases 
exponentially  with  respect  to  the  distance  between  the  two 
points . 

This  section  considers  the  one-dimensional  case.  The 
next  section  will  be  an  extension  of  these  processes  to  the 
two  dimensional  fields.  ^ 
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I Theorem;  The  process  in  2-45  can  be  generated  by  passing 

j white  noise  through  a simple  filter.  The  spectral  density 

of  the  white  noise  is: 

N = 2 aa^ 

o 

and  the  filter  is; 


See  Fig.  11. 


U(t)*WHlTE  NOISe| 

X,n  * MARKOV  PROCESS^ 

H,  * 

1 J w J a J w 

P (<k)  » N„| 

XX  ol 

2 

Fig.  11:  The  Linear  System,  H(jw),  for 

Generating  First-Order,  One- 
Dimensional  Markov  Process 


Proof ; 


Using  2.36  for  one  dimensional  case: 


spectral  density 
at  the  filter 
output 


XX 


(u)  = 


43 


i 


spectral  density 
at  the  filter 
output 


R^(w)  |H{joj) 


The  filter  in  Fig.  11  can  be  represented  by  a differential 
equation,  which  will  be  associated  with  the  process: 

2-47 

where  u(t)  is  white  noise  with  autocorrelation  function 

2-48 

If  a restriction  is  added,  that  the  probability  density 

tt 

function  of  u (t)  is  Gaussian,  the  process  is  called  Gauss' 

ff 

Markov  process,  and  is  completely  described  by  the  auto- 
correlation function.  In  this  case  x(t)  is  also  Gaussian. 
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Another  term  which  is  used  in  connection  with  the 


Markov  Process  is  the  "correlation  time"  point) . This 
time  is  1/a. 

2 . Markov  Sequence  (Discrete) 

Define : 


^-aT  A 
e = p 


T = T’n  n=0,l,2, 


2-49 

2-50 


The  autocorrelation  will  exist  only  for  discrete  points: 


R (n)  = • p'"l  n = 0,1,2,  . . . 

XX  ' 


2-51 


and  taking  the  z transform  of  2-51: 


{1  - P^) 


(1  - P*  z”^)  • (1  - P*  z)| 


2-52 


Theorem:  The  discrete  process  of  Eq.  2-51  can  be  generated 

by  passing  white  noise  through  a discrete  filter.  The 

2 2 

"discrete  spectral  density"  P(z)  of  the  noise  is  a (1  - p ) . 
Fig.  12  shows  the  filter: 


W(|j)  • WHITE  NOISE 

H . — ! 

X(l<)  ^ 

•w 

i-,z- 

j 

R,,  (n)  « <r 

rig.  12:  Discrete  Filter  — The  First  Order  Markov  Sequence 
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Proof : 


Using  2-41  for  one  dimensional  case: 


2 transform  of! 
the  filter 
output 


Rw^(z)-H(z)-H(z"^) 


= a 


(1  - 


1 - pz 


PZ 


L 


transform  of  the  given 
Markov  Process 


Q.E.D. 


The  difference  equation  that  describes  the  filter  is: 


x(k+l)  = p-x(k)  + w(k+l) 

It  is  convenient  to  define 

u(k)  ^ w(k+l) 

u (k)  is  also  white  noise  with  the  same  statistics  as  w(k). 
From  the  last  definition  it  follows: 
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1 


2-53 


(2  2, 
o • (1  - p ) 

if 

n = 0 

R (n  ) = < 

uu  1 

1 0 

if 

n ^ 0 

x(k+l)  = p*x(k)  + u (k) 


2-54 


The  meaning  of  equation  2-54  is  that  the  value  x(k)  depends 
only  on  the  value  of  the  last  step  in  the  past,  x(k-l) . 

3 . The  "Band-Limited",  Continuous  Case 

Another  typical  autocorrelation  function  is: 


R (t)  = 

XX 

^2  i T 1 

a • e ' ' 

• cos  (W  T) 
o 

The  spectral  densi 

ty  will  be: 

R (w)  = 

XX 

2aa^ 

“ 2 1 

(ii)  - u)^)  + a 

This  is  the  case  of  a "band  limited"  signal.  The  process 
x(t)  is  limited  in  a frequency  band  that  is  seen  in  fig.  13. 
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Fig.  13;  "Band-Limited"  Markov  Process 

In  this  thesis,  the  interest  is  mainly  in  discrete  cases. 

The  continuous  case  is  given  here  because  one  can  easily  see 
in  the  continuous  case  the  frequency  band  of  the  random 
signal  and  to  understand  that  the  autocorrelation  function 
in  2-55  represents  an  important  class  of  random  signals. 

In  order  to  make  a complete  summary,  the  dynamic  equations 
or  this  case  are: 


„2  4 2 ^ 2 

6 = a + 0) 

o 


R (t)  = 2aa^6(T) 

uu 


2-58 


2-59 


Now,  the  discrete  case  will  be  discussed  in  detail, 
4 • The  "Band-Limited"  Discrete  Case 

By  using  the  definitions  of  2-49,  2-50: 


- T A 
e = p 


r = T • n 


n = 0,1,2,.,. 


and  define: 


0)  • T = 9 

o 


2-60 


Equations  2-55  has  the  form: 


R (n) 

XX 


= a 


• cos  (9-n) 


n = 0,1,2, 


J 


2-61 


and  the  Z transform  of  2-61: 


49 


P (z)  = (n)} 

XX  XX 


(1-p^)  • [-2-P  cos  9 -t-  (1  + p^)  - Z'^  » cos  9]g^ 

(1  - 2pZ  cos  9 + p^Z^)  (1  - 2pZ~^  cos  9 + p^z 


In  order  to  find  a set  of  different  equations  that  will 
describe  the  given  autocorrelation  function  one  has  to  use 
the  procedure  that  is  described  in  section  D. 

The  first  step  is  to  find  a factored  expression  for 


Pxj^(z)  = Pj^(z)  • Pj^(z”^)  . 

Therefore  P (z)  is  assumed  to  have  the  form: 

XX 


P 

XX 


(z) 


2 2 
a^(l  - p^)  ( 


az  + b 

ZT  2 -2 ‘ 

1 - 2pz  cos  9 + p z 


# 


a*  z + b 

1 - 2pz  cos  9 + 


P z 


2-62a 


The  comparison  of  2-62  to  2-62a  leads  to  a pair  of  algebraic 
equations : 


2 

a 


+ b 


1 + p' 


a*b  = -p  cos  9 


2-63 
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and  the  solution  for  a,b: 


where ; 


6 = 

1 - 

p • cos 

9 

+ 

2 

P 

e = 

1 + 

p • cos 

9 

+ 

2 

P 

2-64 


2-65 


Next,  by  using  2-41,  one  can  see  that  the  random 
signal  discussed  here  can  be  generated  by  passing  white  noise 
through  a discrete  filter: 


WHITE  NOISE 


a Z ' <•  b 


1-2  ^ 1 ^ Z 


Z -,-2 


X (r; 


X hos  the  given  auto  correlation 
(n) 

function  of  e c 2 61 


Fig.  14:  Filter  for  One-Dimensional  "Band  Pass"  Process 


where 


(l  - p^)cr^ 


"ww*'” 


if  n = 0 
if  n 0 
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It  is  convenient  to  define: 


W(z)  = Z”^  U(z)  2-66a 

u(k)  is  also  white  noise  with  the  same  autocorrelation 
function. 


r 

^(1  - pV^ 

if 

o 

II 

c 

II 

1 

1 

0 

if 

n 0 

and  the  filter  will  be: 


X(z) 

wT^ 


a • Z~^  + b 

1 - 2pZ"^cos0  + p^z"^ 


r-6'6^ 


define 


Xj^(z)  = 

X^iz)  = 
From  the  last 

Xj_(k)  = 

x^(k) 


W(z) 

2 -2 

1-Z  2pcos6  + p Z 
definition 

Xj^(k-l)  + 2pcos  0tX2(k-l)  + u(k-l) 


t 


In  matrix  form: 


1 X ■■  (k) 
1 X 

0 -p^ 

'x^(k-l) 

0 

= 

4 

+ 

X2  (k) 

1 2p  cos  0 

1 

^2  (k-1) 

1 

• U(k-l)  : 


and: 


X(z)  = (z)  (a- 2~^  + b) 


x(k)  = a-^2 (k-1)  + b-  X2 (k) 


-2 

= -0  *a‘Xj^(k)  + b*X2(k) 


2-68 


An  Approximated  Model 

Using  the  following  procedure  one  can  derive  an  approxi- 
mated model,  which  is  simpler  than  the  previous  model. 

Define : 


(a-Z'^  4-  b) -W(z)  ^ Z’^*U,(2) 


-67 
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From  the  last  definition: 


u (k-1)  = b-W(k)  + a-W(k-l) 

di 

The  autocorrelation  of  u (k)  is 


!G^(1+P^)  (1-P^)  if  n = 0 

'^^-pcos  0)  (1-p^)  if  n = 1 
0 if  n ^ 2 


Because  R.  (1)  ^ 0,  u is  not  white  noise.  Actually,  u 
u,u,  a a 

a CL 

is  a linear  combination  of  white  noises. 

Now,  an  approximation  is  done 


R (1) 
^a^a 


0 


In  this  case: 


and  now  u (n)  is  white  noise. 

di 

Figure  15  shows  the  correct  and  the  approximated 

noise. 

I 

I 

J 


2-69 


CORRECT 

AUTOCORRELATION 


APPROXIMATED 

AUTOCORRELATION 


Fig.  15:  The  Autocorrelation  of  ^ (n) 

^ a a 


The  filter  that  generates  the  "band  limited"  random 
signal  is  (compare  with  2-66) 

z"^-u  (z) 

X(z)  = 2 T2 

1 - 2pZ  cos  9 + p Z 


Define : 

x^(k)  = -p^-x^Oc-l) 

x^OO  ^ x()c) 


x^  (k) 

0 

-p2 

Xj^(k-l) 

+ 

“ 

0 

. u^fk-l) 

X,  (k) 

1 

2p  cos  9 

X,  (k-1) 

1 

>«(k)  = [o  l] 


X ^ (k) 


X 2 (k) 


This  model  is  simpler  than  the  one  previously  described. 

5 . Second  Order  Markov  Process 

In  the  continuous  case,  the  process  is  defined  by 


define : 


R^^(t)  = • e"“i  ^ I (1  + a|x  I ) 


-aT  A 
e = p 


T = n • T 


aT  = inp  = 9 


It  leads  to  the  discrete  autocorrelation  function; 


R^^(n)  * o^pl"' (1  + 9|nl) 


The  z transform  of  this  process  is 


’’xx'"’  ' 2{R„(n)}  - Pi(x)  + 


where ; 


O0p[4p-(Z+Z  )(l+p)] 


(1  - pZ  ‘)-(l  - pZ)- 


P2(z) 


2 2 


(1  - pZ  (1  - PZ) 


This  process  can  be  generated  by  a combination  of  two  filters 
that  are  forced  by  two  uncorrelated  white  noises: 


H,  (z) 


Xgdc) 


»2  ) 


= X (k)  + X^  (k) 
(k ) a b 


Xb(l.) 


Fig.  16:  Filter  to  Generate  Second 

Order  Markov  Process 


x^(k)  and  Xjj(k)  are  uncorrelated  (because  Uj^(k)  U2^^) 

uncorrelated).  Therefore  the  spectrum  of  x(k)  is  the  sum 

of  the  spectrums  of  x- (k)  and 

ci  D 

Now:  Pj^(z)  can  be  written: 


Pj^(z)  = 


-g^ep  (aZ~^  - b)  (aZ  - b) 


(1  - pZ’-")^(l  - pZ) 


where i 


+ b^  = 4p| 


1 1 

a = + 6^)  e = 4p  + 1 + p^ 


1 1 

17  7 2 

b - ^(e  - 6^)  6 = 4p  - 1 - p^ 


the  filter  H^(z)  is: 


rj  _ (az‘^  - b) 

^ (1  - PZ 


(z) 

u^TzT 


1 f 

2 

1 1 

-9pa^ 

if 

3 

II 

0 

1 

L 0 

if 

m ^ 0 

the  state  variable  expression  that  follows  Hj^(z): 


X,  (K) 


1 r 

-p2  I x^(k+l) 


Xj  (k-1) 


(k)  = (-P  a -b) 


xj_(k) 


X2(k) 


and  for  H2(z): 


H2 (z)  = 


1 

1 - pZ" 


Xb(z) 

wJTzT 


define : 


U2(k)  » w^Ck+l) 
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6.  Conclusions 

1)  It  was  shown  in  this  section  that  some  stochastic 
processes  can  be  generated  by  passing  white  noise  (or  a . 


combination  of  white  noises)  through  a linear  filter  of 
order  P. 


2)  Also  it  is  seen  that  the  order  P is  a finite 
number  [in  the  first  two  examples  P was  1,  in  the  last 
two  examples  P was  2] . 

3)  Pj^,  the  order  of  the  process,  tells  the  number 
of  "neighbors"  near  the  point  k,  upon  which  the  value  X(k) 
depends . 


F.  STATE  VARIABLE  MODELING  OF  TWO  DIMENSIONAL  FIELDS  BY 
"FILTER  RESPONSE  METHOD" 


1.  Introduction 

The  starting  point  of  the  modeling  procedure  will 
be  an  extension  of  the  Markov  Process  to  a two  dimensional 
case. 

It  is  emphasized  that  the  technique  used  in  this 
section  is  good  only  for  separable  autocorrelation  functions. 
In  Section  H another  modeling  method  which  is  valid  for 
any  homogeneous  random  field  will  be  discussed.  [A  compari- 
son between  the  two  methods  is  given  in  Section  I] . 

2 . Model  For  First  Order,  Continuous,  Markov  Process 
Given  a two  dimensional  field,  with  an  autocorrelation 

function ; 


R (t 

XX 


1/T2] 


2 ’“l 
a e 


2-73 
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I 


is  the  correlation  factor  in  the  t^  direction 
a.2  is  the  correlation  factor  in  the  t2  direction 


Fig.  17:  System  Coordinates  For 

The  Field  in  2-73 


The  power  spectral  density  function  (using  2-46,  2-2 


(2a,  ) (2a 


(o)^  + a^^  ) (o)2  + ct2  ) 


This  autocorrelation  function  might  be  generated  by  passing 
white  noise  through  a two  dimensional  filter: 


u ■ wMiTe  Noise  h . ! 

r,y(«,_Uj)  • 40|  


<"|  "j  > • OlVEN  I 

[of  Ey 


Fig.  18:  Generation  of  First  Order,  Two 

Dimensional,  Markov  Process 


The  proof  is  by  using  2-36.  One  can  see  that  the  spectral 
density  at  the  filter-output  is  the  same  as  in  2-74. 
Therefore: 


U (u)^,0J2) 


{ j(jj^  + ( ju)2  + ^2^ 


Define : 


M(a)^,o)2) 


U (a)^,u)2) 


Then: 


dtl 


-U2M(t^,t2)  + '^(tj_,t2> 


N (ojj^,(J2) 


M(a)j^,W2) 


dN(tj^,t2) 

3t"i 


-aj^N(tj^,t2)  + M(tj^,t2> 


and  in  matrix  form: 


62 


ft 

(-■ 

= n 

1 = 

0,1,2,  . . . 

^2 

= kT 

k = 

0,1,2,  ... 

'^1 

= nT 

n = 

^2 

= mT 

m = 

0/1^2/  *** 

has 

now  the 

form: 

R 

XX 

(n,m)  = 

n 1 |m  1 

°2 

n = 0,1,2,  ... 
in  0,1,2,  ... 

2-78 


and  the  two  dimensional  2 transformation  of  2-78: 


o^a  - (1  - P2^) 


(1  - (1  - (1  - p2^2  ) (1  * ^2^2^ 


2-79 


This  discrete  autocorrelation  function  can  be  generated  by 
passing  white  noise  through  a two-dimensional  discrete 
filter: 


64 


Fig.  19:  The  Discrete  Area 
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It  is  convenient  to  define: 


W(2^,22^  - 22’^  U(Zl,Z2^ 

and  the  filter  has  not  the  form 


can  continue  in  two  directions : 
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Direction  1;  (Model  1) 


From  2-80: 


^(Z2»Z2)'[1  — — z 2 0 2 2^  z 2 ^1^2^  ~ ^ 


and  a difference  equation  can  be  written; 


x(k+l,i+l)  = pJ_■xOc+l,^)-p^02•x(k,^)  +P2-x(k,i+l)  +x(k,Jl) 


R^(n,m)  = 


2 2 

(1  - Pj^  ) (1  - p2*”)  if  and  m=0 


if  n?^0  Of  m?^0 


6 


Direction  2:  (Model  2) 


Define 


N(z^,Z2) 


1 - Z,  -P, 


M(z^,Z2)  = X(z^,Z2) 


From  2-83 


N(z^,Z2)  = (z^,  Z2)  + Z^”^  U(z^,Z2) 


N(k,Ji)  = p^N^(k,£-l)  + u(k,£-l) 

Note:  from  Eq.  2-78  it  is  obvious  that  Z2  corresponds  to 

the  k direction  and  z^  corresponds  to  the  I direction. 
Substituting  2-83,  2-84  into  2-80: 


z_”^N(z, ,z,; 
1 - Z2  P2 


M(z^,z^)  = Z2~^P2M(Zj^,Z2)  + Z2~^N(Zj^,Z2) 


M(k,£)  = P2M(k-l,£)  + N(k-1,£) 
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By  changing  indices  of  2-85,  2-86  one  can  write  these  two 
equations  in  a state  vector  form 


N(k,Ji+l) 


M(k,il) 


N(k,£) 


W(k,)l) 


x(k,Jl)  = ( 1 0 ) • 

M(k,il) 

N(k,Jl) 

By  comparison  2-87,  2-88,  to  2-81 


M{k,L)  = '<(’^,1) 


N(k,L)  = X(k+1,«,)  - P2-X(k,£) 


It  is  obvious  that  model  1 and  model  2 describe  the  same 
field. 

An  interesting  property  of  this  model  : 


E{M(k,il)N(i,  j)  } = 0 

for  any  k,2.,i,j. 
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Proof : 


E{M(k,Jl)N(i,  j)  } = E{x  (k,Jl)  [x  (i+1,  j)  - P2  x(i,  j)  ] } 

= E{x  (k,&)  [p2  x(i,  j)  +u(i,j)  -P2x(i,j)]} 


u(i,j)  is  white  noise  and  therefore  it  is  uncorrelated  to 
X (k , £ ) . Therefore : 

E{M(k,£)N{i, j) } = 0. 

Q.E.D. 


4 . Model  For  "Band  Limited’*,  Discrete  Case 

Equivalently  to  the  previous  discussion,  the  two- 
dimensional,  "band  limited."  discrete  Markov  Process  is 
defined  by  the  autocorrelation  function: 


•y  |t%|  )ml 

R^^(n,m)  = [p3_  ' cos  (9  ^^n)  ] [p2  ' 'cos(92m)J 

m=  0,1,2,  ... 


2-91 

The  z transform  of  this  autocorrelation  function: 


P 

xx 


Z2) 


a^(l  - Q^)  (1  - P2^) 


B (z~, Z2) 


2-92 
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L 


where 


A(z^,z 


B(Zj^,z 


A(z^,Z2 
A(z,  ,z. 

X 4 

where ; 


2-1  2-1 
,)  = [-Zj^cos9^+(l+pj^  )-z^  cos6^]  [-Z2OOs02+(1+O2 


,)  = [ (l-2pj^Zj^cos0j^+O^^z^^)  (l-2p^z^~^cos0^+p^z^~^)  ] 


[ (l-2p2Z2Cos02+p2^Z2^^  (l-2p2Z2~^cos02+P2^Z2  ^)J 


) can  be  written  in  a form; 


2-93 


2-94 


,)  = [ (a^z^"^+bj^)  (a2Z2+bj^)  ] [ (a2Z2~^+b2)  (a2Z2-HD2)  ] 2-95 


2 2 2 
+ b^^  = 1 + p^^ 


2-96a 


2 2 2 
a2  + ^2  = ^ ^ ^2 


2-96b 


aibi 


= -pj^cos9j^ 


2-96c 


*2^2 


-p2COS02 


2-96d 
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> 


The  solution  for  the  a's  and  b’s: 


Now:  This  random  process  can  be  generated  by  passing  white 

noise  through  a discrete  filter  of  the  form: 


2 


Fig.  22:  Filter  to  generate  "Band  Pass"  yandom  field. 


where : 


jc3(l-p^  )(l-p_)  ifn=0  and  m = 0 

R (n,m)  = ( ^ > 

I 0 if  n ^ 0 or  m ^ 0 


It  is  convenient  to  define: 


Vt(z^,z^)  = ^l”^‘^2~^  ■ 


u(k,)l)  is  also  white  noise  with  the  same  statistics  as  w(k,i-) 


the  filter  has  the  form: 


X(Zj_,Z2) 


^ (^2^2  ^ * ^2^  ^1  ^^2  ^*U(Zj^,Z2) 


(1  - 2Pj^z^  cos  + z^  )(1-2p2Z2  cos  02  + P Z2  ) 


Now  the  following  definitions  are  done: 


N^(Zj^,Z2)  = "Pi^*Zi~^‘N2  (Zj^,Z2) 


^2^^1'^2^  = ^ ^ ^ 

1 - 2Pj_z^  -^cos  9j_  + Pj_^Zj_ 


^1^^1'^2^  = ^ + bj^)  W2  (z^,  Z2; 


From  the  last  three  definitions  one  can  write  a set  of 
difference  equations: 


(k,£) 


N2(k,lt) 


0 -P, 


(k,il-l) 


1 2pj^cos  0J  |N2{k,K.-l)J 


2-102 
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Xj^(k,8.)  » -pj^  *aj^N^(k,i)  + bj^N2(k,i) 


Now,  other  definitions  are  done: 


2 -1 

= -p2  *22  •M2(zi,Z2) 


Zm  *^1  ^Z-i/Z^) 

= - --  - 2 : 
1 - 2/2*Z2  -cos®  2 ■^'^2  ’^2 


From  these  definitions  it  follows: 


X(Zj^,Z2)  = ^®2’^2  ■*■  ^2^  ‘ ^2^^1'^2^ 


Mj_(k,ii) 


M2(k,£) 


^ 2P2COS  02  M2(k-1,£) 


(k-1,  £) 


x(k,£)  * -P2"  *a2*Mj^(k,£)  + b2*M2(k,£) 
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Equations  2-102,  2-103,  2-105  can  be  written  together: 


1 


Using  2-103 


Approximated  Model 


One  can  derive  for  the  "band  limited"  case  an  approximated 
model,  which  has  a simpler  structure.  The  starting  point  is, 
again,  2-101* 


Define ; 


Ua(Zi,Z2^  = (a^z^"^  + (a2Z2~^  + b2>  u(z^,Z2)  2-108 


From  eg.  (2-87)  one  can  see  that  u (k,£)  is  the  noise- 
input  to  the  filter  that  creates  the  "band  limited"  random 
field.  Now,  what  are  the  statistics  of  u (k,)i)?  Equation 
(2-108)  tells  us  that  u (k,J,)  is  a linear  combination  of 

di 

white  noises; 


u^(k,i)  = a^^a2U(k-l,Jl-l)  + b^^a2U(k-l,  + aj^b2u(k,Jl-l)  + bj^b2u(k,£) 

2-109 


Using  eq.  (2-108),  together  with  the  definitions  of  the  a's 
and  b's  in  eq.  (2-96),  the  autocorrelation  of  u (k,Jl)  can  be 

3 


u,u  (n,m) 

9 d 


(1+Pj_2)  (1+P2^)  (1-Pi^)  (1-P2^) 

(a  nunfcer  that  is  not  zero) 


if  ra  = 0 and  n = 0 


if  m < 1 and  n < 1 


if  m > 1 or  n > 1 

2-110 
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For  example,  the  calculation  for  R (0,0)  is  shown: 

R „ (0,0)  = E{u^^0c,Jl)} 

Va  ^ 

= [(a^a2)^  + (bj^a2)^  + + (bj^b2)^]  • [ (1-p^^^)  (1-P2^)  ] 

= (a^2  + b^^)  (a2^  + b2^)  (l-p^^^)  (1-P2^) 

= (1-p^^)  (1+P2^)  (1-p^^)  (1-P2^)  2^111 


Here  is  the  place  to  make  an  approximation:  One  can  see  that 

u white  noise,  because  R ^ (n,m)  0 when  n or  m 

a a a 

are  equal  to  one.  Now,  arbitrarily,  it  is  decided  to  let 

this  value  ^ zero  (for  n or  m equal  one^.  In  that  case  u is 

3i 

white  noise: 


«u  u 
a a 


(1+Pj^2)  (1+P2^)  (1-p^^)  (1-P2^) 


if  n = 0 and  m = 0 

> 

if  n 0 or  m ^ 0 


2-112 


The  filter  to  create  the  random  field  will  be: 


X(Zj^,Z2) 


. ... 

TT  2 — IT — 5 — IT  2-113 

(l-2p^z^  COS  ^1  ) (l~2p2Z2  "^cos  ^2 '^2^*2  ^ 
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Comparing  2-113  with  2-101,  one  can  write  for  2-113  a state 
variable  model  that  will  be  similar  to  eq.  2-106: 

Define : 


From  these  definitions  it  follows: 


2-114 


8 


By  comparing  the  difference  equation  that  follows  from  2-113 
to  the  structure  of  2-114,  2-115,  it  is  found: 


M2  ik,i) 

= x(k,^) 

N2  (k, «.) 

= x(k,Jl+l)  - 2p2"Cos  02' 

x(k,Jl)  + P2*x(k,£-1) 

Mj_(k,i) 

= -P2^-x(k-l,il) 

(k,£) 

2 

= -p^  •x(k,2.)  - 2p2‘COS 

02-x(k,2.-l)  + P2*x(k,2.-2) 

5 . Model  for  Second  Order  Markov  Process 

The  extension  of  the  one  dimensional  case  to  the  two 
dimensional  case  is: 


R (n,m) 

XX 


= Pi 


n 


(1  + 0. 1 n I ) p. 


m 


(1  + 


92|m| ) 


and 

the 


following  very  similar  procedures  to  all  previous  cases, 
state  space  model  in  this  case  has  the  form: 


I 
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4p^  + 1 + 


1 1 
111 
T (^1  + ) 


4P2  + 1 + P2- 


1 1 

1 '2'  , ^ Iv 

2 ^^2  ^ *^2  ^ 


4P,  - 1 - P,' 


1 1 
I ■ '^1  > 


4P2  - 1 - P2' 


1 1 
^ fp  ^ ^2. 

2 ^^2  ■ *^2  ^ 


u^u^(m,n) 


if  m = n = 0 


if  m 7^  0 or  n 7^  0 


R.  „ (m,n) 


^2^^2 


(1-pj^^)  (1-P2^)a^  if  m = n = 0 


if  n 0 or  m ^ 0 


G.  THE  STATE  SPACE  STRUCTURE 


1.  Introduction 


Section  F showed  that  one  can  arrive  at  a Discrete- 


Space-Model  for  linear  Image  Processing,  that  has  the  form: 


. u(k,«,) 


M(k+l,il) 

2 

^1  : ^2 

1 

1 

• * T - - ~ 

1 

0 

M(k,)l) 

+ 

®1 

N(k,  Jl+1 

A3  : A^ 

N(k,£) 

^2 

X (k,  Jl) 


M(k,Jl) 

N(k,Z) 


where : 

k:  An  integer  valued  vertical  coordinates 

Z:  An  integer  valued  horizontal  coordinates 

M:  A vector  which  conveys  information  vertically 

N:  A vector  which  conveys  information  horizontally 

u:  A vector  that  acts  as  an  input 

x:  A vector  that  acts  as  an  output. 

Aj^,  A2/  A^/  B2/  C2  are  matrices  of  appropriate 

dimensions.  Boundary  conditions  N(k,0),  M(0,JI)  are  also 
inputs.  u(k,)l)  is  specified  externally. 

This  section  will  summarize  the  properties  of  the  State- 
Space  model.  The  discussion  is  based  on  Ref.  [12] . 
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■ .ly  iLiiii 


T 


2.  Realization  of  a Discrete  Filter  By  The  State 
Space  Structure 

The  following  realization  result  has  been  demonstrated; 
Given  an  arbitrary  two  dimensional  digital  filter  whose  2-D 
z transfer  function  is: 


n(z^,Z2) 


_ m _ n , _ _ m- 1 _ n , 

Z,  *Z_  + a,^*Z,  ’Z_  + ...  a 

1 2 10  1 2 mn 


2-117 


then  matrices  A^i  A^  and  vectors  exist  with 

dimensions  A^  - mxn,  A2  — mxn,  A^  — nxm,  A^  — nxn,  — mxl, 
B2  - nxl  such  that  Eq.  2-117  may  be  put  into  the  form  of 
Eq.  2-116. 


In  addition,  the  following  canonical  form  for  Eq. 
2-116  has  been  developed:  If  the  denominator  d(z^^,Z2)  of 

2-117  factors  as  follows  (such  a factorization  is 
denoted  "doubly  factorable')): 

d(z^,Z2)  = d^^(Zj^) -d2^(z2)  + d^^(Zj^)  *d2^(Z2) 


then  the  canonical  form  for  the  A^  (i  = 1,...4)  of  Eq. 
2-116  is: 
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a 


where  the  coefficients  a.,  b. , c.  and  d.  are  determined 

111  1 

from  d(z^,Z2)  by 


/ m™  i , _ m™  2 * \ ^ « n ™ 1 , 

* (b  z . + b 1^1  +.*.b^  + . • • c 

m 1 m-  i 1 1 " i 


2.  The  State  Transition  Matrix 


Definitions: 


A ^ 


^ A2 


A3  A^ 


2- 


B ^ 


B, 


2- 


C 5 ic^  C2 


2- 


M(k,£) 


T(k,£)  = 


2- 


N(k,£)| 


87 


118a 


118b 


118c 


118d 


Then 


Now, 


T' (k, £)  = 


M(k+1,£) 


N(k,£+l)l 


2-118e 


T'  (k,£)  = AT(k,£)  + Bu(k,£) 


2-119 


y(k,£)  = CT(k,£) 


2-120 


for  k > 0,  £ > 0,  the  next  definitions  are  done: 


,1.0 


2-121a 


,0,1 


Ai  A2 


2-121b 


^k,  £ _ ^1,0  ^k— 1,£  ^ ^0,1  j^kfZ—1 


2-121C 


A»'°  =I 


,-k,£ 


. 0 
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A(k,£.)  is  called  "State  Transition  Matrix".  Examination 
of  these  definitions  leads  to  formulas  of  system-response 
that  are  similar  to  the  one-dimensional  discrete  case. 

Table  1 compares  the  formulas  for  the  one 
dimensional  case  to  the  two  dimensional  case. 


3 . Properties  Of  A 


k,  I 


1,  A = 


A,  A- 

A,  A., 

1 2 

1 2 

= 

+ 

A,  A. 

0 0 

3 4 

0 0 


A3  A4 


A = A^'°  + A°'^ 


2-122 


2. 


Ak,0 


a1,0  ^k-1,0  ^ ^0,1  ^k,-l 


= A^.O  ^k-1,0 

Thus , 
and 


aJ^,0  ^ 

Z 

} 

2-123 


2-124 


3.  I = 


I 


0 


0 


I 


89 


Table  1;  Comparison  Between  the  Response  of  Two 
Dimensional  Systems  to  One  Dimensional 
Systems 
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where  I is  the  identity  matrix  with  appropriate  dimensions 
Thus , 


or : 

2-125a 
2-125b 


e:  T®'  ^ 0 


0 


Briefly, 


jO'! 


= 0 


2-126 


In  the  next  discussion  we  shall  use  some  figures 
for  better  understanding.  The  figures  include  arrows  that 
are  going  from  one  pixel  to  its  neighbor.  For  excimple: 

The  arrows  in  Fig.  23  means  that  the  values  N(k,i) 

contribute  the  values  N(k,i+1)  and  M(k+l,i). 


Fig.  23:  The  Propagation  of  the  States  M and  N. 

[The  values  M(k,i),  N(k,i)  contribute 
the  values  N(k,i+1)  and  M(k+l,i)]. 


4 . General  Response  Formula 

Lemma:  Let  u(k,i)  be  zero  for  all  (k,i). 

M(0,i)  « N(k,0)  = 0 for  (i,j)  0 . 
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(only  M(0,0),  N(0,0)  are  not  zero). 


Then: 


T(k,Jl)  = T(0,0) 


2-127 


Proof:  (By  induction) 

T(0v0)=  I-T(0,0)  = T(0,0) 

Now  assume  it  is  correct  for  any  (0,0)  £ (i,j)  £ (k,?,)  , then; 

I M(k,^)| 


T(k,)l)  = 


N(k,Jl)| 


A^  M(k-1,«,)  + A2  N(k-l,^)  + B^-0 


A^  M(k,£-1)  + A^  N(k,£-1)  + B^-O 


Ai  A2 

0 0 

ss 

T(k-l,JO  + 

0 0 

A3  A^ 

= A^'°  T(0,0)  + A°'^  T(0,i 


a’^'*'  T(0,0) 


Q.E.D. 


Figure  24  shows  the  propagation  of  the  field  due  to 
M(0,0)  and  N(0,0) . 
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Fig.  24:  The  propagation  of  states  due  to  M(0,0)  N(0,0) 


Equation  2-127  will  be  used  for  the  next  two  cases: 

Effect  of  M (0 , j ) : 

Assume:  M(0,j)  is  the  only  non-zero  boundary  condition 

and  that  all  inputs  are  zero. 


Using  Eq.  (2-118d) : 
T(0,j)  - 


M(0,  j) 
0 


9 


Therefore  T(0,j)  can  be  used  as  initial  condition  for 
2-128  : 


2-128 


Effect  of  N (i, 0)  : 

Similarly  to  M(0,j)  the  effect  of  N(i,0)  is: 


Effect  of  u (k,  ^)  : 

Assume:  u(i,j)  for  some  (i,j)  < (k,il)  is  the  only 

non-zero  input.  All  boundary  conditions  are 
zero. 

Then: 


T(i+1, j) 


M(i+1, j) 

N(i+1, j) 

Aj^M(i,j)  + A2N(i,j)  + B^u(i,j) 

A2M(i+l,j-l)  + A2N(i+l,j-l)  + B2U(i+l,j-l) 


T(i+1, j) 


A^'O  + A2*0  + B^-u(i,j) 
0 


and: 


u(i,  j) 


I T(i,j+1) 


0 


B, 


u(i, j) 


T(i+l,j),  T(i,j+1)  might  be  substituted  in  2-127  as  boundary 
conditions.  Therefore,  by  using  superposition: 


Fig.  25  shows  the  effect  of  u(i,j). 


I 
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Proof ; By  superposition  of  the  effects  of  all  inputs  and 


boundary  conditions . 


Q.E.D. 


5 . Characteristic  Function 

In  the  one  dimensional  case  the  eigen  values  of 
the  system 


x(k+l)  = Ax(k) 

are  defined  as  those  values  which  satisfy  the  algebraic 
equation; 

Ax  = Xx 

or  the  characteristic  equation 
|XI  - a|  = 0 

Now,  for  the  "State  Space  Structure"  let's  define  in  the 
seune  way:  given: 

M(k,Jl) 

N(k,)l) 

i 

1 


M(k+l,H) 

^1 

^2 

N(k,A+l) 

^3 

^4 

define  operators  E,  F such  that: 


^1 

^2 

M(k,£) 

EM(k,i) 

^3 

^4 

N(k,i) 

FN(k,il) 

or: 


E -I-A^ 

^2 

M(k,£) 

^3 

F*I-A^ 

• 

N(k,£) 

0 2-131 


To  have  a nontrivial  solution  of  this  equation  we  require 
that  the  matrix  in  2-131  is  singular.  Therefore  the 
determinant  should  be  zero. 


2-13l§ 


6.  Stability 

The  stability  criteria  of  Huang  [24]  can  be  generalized 
in  a straightforward  manner  to  systems  represented  in  state 
variable  form.  This  generalization  allows  the  use  of 
standard  one-dimensional  routines  in  the  determination  of 
two  dimensional  system  stability. 
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The  stability  criteria  will  be: 

1)  The  eigenvalues  of  the  matrices  A^,  are 
less  than  one  in  magnitude. 

2)  The  eigenvalues  of  the  complex  matrix: 

A4  + A3(z^  2-132 

are  all  less  than  1 in  magnitude  as  varies 
about  the  unit  circle  (|z^|  =1).  If  any 
eigenvalues  of  Eq.  2-132  has  magnitude  greater 
than  one  for  = 1,  then  the  system  is 

unstable. 

H.  MODELING  BY  USING  OPTIMAL  ESTIMATION  THEORY 

So  far  we  have  seen  a technique  for  modeling  random 
fields  that  uses  z transformation,  properties  of  linear 
filters  and  special  types  of  correlation  functions  (separable) . 
That  method  was  called  "Filter  Response  Method".  In  this 
section  optimal  estimation  theory  is  used  to  solve  the  modeling 
problem.  The  advantage  of  this  method  is  that  it  is  not 
limited.  For  comparison  between  the  two  methods  (see  Section 
I)  . 

1.  The  Basic  Principle 

Suppose  a discrete  random  field  x(k,£)  is  given  and 
assumed  to  be  homogeneous  (stationary) . Given  the  values 
x(k,l)  at  all  points  n(k,il),  the  problem  is  to  estimate  the 
value  x(k,2,).  In  other  words,  if  k(k,4)  is  the  optimum 
estimate  for  x(k,Jl),  then: 


100 


is  minimized.  x(k,)l)  will  be  the  "linear  least  square 
estimate"  of  x(k,)l). 

Substituting  2-133  in  2-134,  differentiating  with 
respect  to  each  ^ , and  setting  each  derivative  equal  to 
zero,  we  obtain  the  following  set  of  simultaneous  equations 
for  the  unknowns  a . . : 

E{[x(k,Z)  - x(k,H)]  x(i,j)}  = 0 for  all 

(i,  j)  tO. 

2-135 


which  says  that  the  coefficients  a.  . must  be  such  that 

1/3 

the  estimation  error  x(k,il)  - x(k,Jl)  is  statistically 
orthogonal  to  each  x(i,j)  that  is  used  to  form  the  linear 
estimate.  This  is  known  as  the  orthogonality  principle,  in 
linear  least  square  estimation. 

Let  D represent  the  following  collection  of  pairs 

(i/ j) : 


D = {(0,1),  (1,1),  (1,0)} 


2-136 


Now  comes  the  definition  of  a first  order  process. 
Definition  1; 

A random  field  will  be  called  first  order  Markov  if 

the  coefficients  a.  . in  2-133  are  such  that  x is  of  the 

1/3 

form: 
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t-l  i 


►X 


Fig.  27:  Definition  of  Reg  D 


x(k,£)  = S a.  . x(k-i,2.-j)  2-137 

(i,j)CD(k,£) 

(all  other  a's  are  zero). 


That  is,  the  least  square  estimate  of  x{k,£)  in  terms  of 
n(k,Jl)  is  the  same  as  that  in  terms  of  only  the  three 
immediate  neighbors  left  and  above  the  point  (k,)l). 

Substituting  2-137  in  2-135,  the  folxowing  conditions 
for  the  Markov  field  must  be  satisfied; 


E{[x(k,)l)  - Z a.  . x(k-i,i-j)]  x(p,q)}  = 0 

(i,j)£D(k,)l) 


2-138 


for  all  (p,q)  < D. 
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B 


For  a Markov  field  the  coefficients  a.  . must  be  such  that 

If  D 

2-138  is  satisfied  for  all  (p,q)  c J2(k,i).  In  particular, 
2-138  must  be  satisfied  for  the  following  values  of  (p,q) : 


(k-l,Jl),  (k-l,ll-l),  (k,H-l). 


Substitutiiig  the  values  of  (p,q)  in  2-138,  the  following 
equations  are  obtained,  for  a. 


a,  „ R (0,0)  + a,  , R (-1,0)  + a.  . R (-1,1)  = R (-1,0) 

1 , 0 XX  ' 1 , 1 XX  0 , 1 XX  XX 


a,  ,,  R (1,0)  + a,  , R (0,0)  +a.  , R „(0,1)  = R (-1,1] 
1,0  XX  ' 1,1  XX  0,1  XX  ' XX 


“1,0  * “1,1  ''xx<°'-l>  * “0,1  \x<°'°>  - "xx*"--!* 


2-139 


where  R (a, 6)  = E{x(k,2,)  X (k+a , H+6 ) } . For  system  coordinates 
in  this  problem  see  Fig.  19. 

Example  1 

Given 


R (n,m) 

XX  ' 


2-140 


by  substituting  2-140  into  2-139, 
a . , is : 


If  3 


the  solution  for  the 
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“o,l  = 


“l,l  “^1  ^2 


“1,0  = ^2 


and  by  substituting  into  2-137: 


2-141 

Note  I The  discussion  in  this  section  will  be  limited  main- 
ly to  "one-side  process"  (causal  difference  equation) 
in  order  to  compare  it  to  the  method  of  modeling  that 
is  described  in  Section  F.  In  the  end  of  this  section 
there  will  be  some  examples  of  non-causal  models. 

2.  The  Modeling  Error 

Definition  2:  The  modeling  error  is  the  difference 

between  the  true  value,  x(k,2.)  and  the  estimate,  x(k,il). 

A 

u(k,il)  = x(k,Jl)  - x(k,H) 


= x(k,£.)  - E x(k-i,?.-j)  a.  . 2-142 

(i,j)en 


It  is  obvious  that: 


x(k,l)  = x(k,£)  + u(k,il) 
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x(k,2.) 


2-143 


Z a.  . x(k-i,H-j)  + u(k,A) 


The  error  u(k,£)  creates  also  a random  field.  The  question 
is  what  kind  of  random  field?  We  are  concerned  in  the 
variance  and  uhe  autocorrelation-function  of  this  error, 
a)  The  Variance  (of  the  modeling  error) 

Using  2-137: 

Q = E{u^(k,i)} 

/s 

= E{  [x(k,^)  - x(k,Z)  } 

= E{x^(k,£)  - 2x(k,H) -x^k.^)  + 

= E{x^(k,H)  - X (k,  ?,)  . X (k,  K,)  + x 

The  zero  in  the  last  term  is  set  by  using  the  orthogonal 
principle  (Eq.  2-138). 

Therefore ; 

Q = E{x^(k,i) 

and  recall  that: 

R(i,j)  = E{x  (k-i,  «,-j ) x(k,Jl)}  = autocorrelation  function 


- [ Z a..  ,.x(k-i,«,-j)  ] x(k,^)} 

(i,j)cfi 
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Q = E(u^(k,^)}  = R(0,0)  - Z a.  . R(i,j) 

(i,j)en 


2-144 


b)  The  Autocorrelation  (of  the  modeling  error) 
Theorem:  The  modeling  error  creates  a random 

field  that  is  white  noise,  e.g. : 


R (n,m)  = 


if  m = 0 and  n = 0 


0 ifm^O  or  n ^ 0 


Proof:  By  using  2-143,  Eg.  2-137  can  be  written  in  a 

non-recursive  way  as  follows: 

k Z 

x(k,i)  = Z Z 6 u(m,n) 

ni,n 


Z' 


z 

I Y 

j = l 


0,  j 


Initial 


x(0,Jl-j)  + 
conditions 


But  using  the  orthogonality  principle; 


E(u(k,)l)  x(i,j)}  * 0 (i,j)  e Q 2-146 


and  especially; 
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0 


2-147 


E{u(k,i)  x(k,^)  } = 

substituting  2-145  into  2-147,  the  result  is  obvious: 
E{u(k,£)  u(m,n)}  = 0 {m,n)  c 

Q.E.D. 

Note:  As  previously  determined,  this  discussion  concerns 

causal  models.  In  this  case  we  proved  that  the 
random  process  is  forced  by  white  noise.  That  is 
not  the  case  in  a non-causal  model. 

Example  2 

Problem:  Find  the  variance  of  the  modeling  error  for  the 

model  of  Excimple  1. 

Solution : 


Q = 1 - R(0,1)  - R(l,l)  - Q R(1,0) 

■ I - ^ - p/ 

Q - (1  - (1  -P2^)  2-148 

Conclusion:  For  the  autocorrelation  function  of  2-140 
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the  model  is : 


x(k+l,)l+l)  = (ktl,  il)  + Pj^p2X  (k,  Z)  +P2x(k,£+1)  +u(k,£) 


E{u^(k,il)  } = Q = (1  - p^^)  (1  - P2^) 


2-149 


Now,  comparing  2-149  to  2-81,  2-82  it  is  seen  that: 

For  the  First  Order  Markov  Field  the  orthogonality 
principle  (Minimum  Mean  Square  Error)  leads  to  the 
same  model  as  the  "filter  response  method". 

3 . Advantages  Of  The  "Orthogonality  Principle"  Method 
There  are  three  distinct  advantages  of  this  method: 
a)  This  method  can  be  applied  to  non-separable  autocorrela- 
tion functions.  This  is  impossible  to  do  with  the 
Linear  Filter  Response  method. 

b.  This  method  can  be  extended  to  non-causal  models.  The 
technique  is  very  similar  to  the  causal-modeling  tech- 
nique. The  only  differences  are  that  the  modeling 
error  is  not  white  any  more  and  that  the  solution 
doesn't  exist  for  all  cases. 

c.  It  is  optimal  in  the  sense  that  the  modeling  error, 
u(k,J.),  is  minimum.  It  was  not  proved  that  the  method 
of  Section  F,  Linear  Filter  Response,  leads  to  an 
optimal  model.  [We  see  it  only  for  one  special  case. 
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the  first  order  Markov  field. ] For  further  discussion 
of  this  point  see  Ref.  23. 

The  next  examples  will  highlight  these  three  advantages. 


Example  3 (for  advantage  a) : 
Given:  R 


XX 


- P 


2 2 
n +m 


2-150 


Find:  A linear  model  for  the  field  that  uses  region 

D (defined  in  Eq.  2-136  and  described  in  Fig. 
27) 


Solution:  As  in  Example  1,  the  set  of  three  algebraic 

equations  of  2-139  has  to  be  used: 


«xx(0' 

0) 

+ 

+ 

I 

‘^1,0 

0) 

+ 

Rxx(0,0) 

+ 

0 

= 

-1) 

+ 

+ 

“o,i 

and  by  substituting  2-150: 


/t 

1 0 n • 

L 1 

X M M 

1,0 

M 

/I 

PIP 

0 

“i,i 

s 

p 

p ^ p 1 

“o,i 

p 
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I The  a’s  are  solutions  of  these  three  equations. 

Note ; The  process  in  Eq.  2-150  is  not  first-order,  if  we 
* use  definition  1.  There  is  not  a finite  set 

of  elements  on  which  the  point  x(k,£)  depends. 
Therefore  the  choice  to  model  the  process  by  region 
D was  arbitrary  in  this  exeunple. 

Example  4:  (for  advantage  b) : 


Fig.  28:  A Non-causal  Region  For 

Modeling  (in  exaunple  4) 


Problem;  Find  the  non-causal  representation  of  x(k,Z) 

by  using  the  dotted  elements  of  Fig.  28,  for  the 
autocorrelation  function 


"xx'"'”’ 


n m 
'P2 
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Solution;  By  using  Eq.  2-138  for 


(i,j)  = {(-2,-2),  (-1,-1),  (-1,1),  (-1,2)} 


we  have  the  next  set  of  equations; 


R^,(0,0) 

^^-2,-2 

R^x(O'O) 

• 

“-1,-1 

R^x(O.O) 

R (-1,1) 

XX 

“-1,1 

R (0,0) 

XX 

“-1,2 

XX 


XX 


XX 


XX 


(-2,2) 


(-1,1) 


(-1,1) 


(-1,2) 


1 P P P ^P  P , ^ 

a 

^ 12  12  1 

-2.-2 

1 2 

PP,  1 P^' 

a , 

P ,P 

12-^  1 12 

-1,-1 

1 2 

p 3 p p ^ 1 Pj_  P2 

X X 

“-i.i 

P^P2 

“’l"  "l^‘=2  "l"2  '■ 

“-1.2 

2 2 
Pi  P2 
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The  solution: 

“-2,-2  ^ “-1,2 
“-1,-1  = “-1,1 

It  is  not  surprising  that  a_2  _2  2 zero. 

The  given  autocorrelation  function  in  this  excimple  is  a 
first-order  Markov  process. 

Example  5; 

This  example  compares  between  the  two  methods  of  modeling 
(by  optimal-estimation  and  by  "filter-response-method"). 
The  results  turn  out  to  be  different.  In  this  example  we 
consider  a second-order  process.  The  optimal  estimation 
approach  for  the  two-dimensional  case  requires  solving 
eight  algebraic  equations,  which  is  complicated. 

Therefore. this  example  deals  with  a one-dimensional  auto- 
correlation function.  The  extension  to  a similar  two- 
dimensional  problem  is  straightforward. 

Given:  The  one  dimensional  "Band  Pass"  autocorrelation 

function  (see  Chpt.  II,  Section  E,  parts  3,4). 

R (n)  * cos  (9n) 

xx 


= 0 


P1P2 

1 + p. 
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Find ; A model  that  describes  this  process  by: 

a)  optimal  estimation  approach. 

b)  by  filter  response  method. 


Solution: 

a)  In  the  optimal-estimation-approach  one  has  to 
determine  the  order  of  the  model  (order  of  the  difference 
equation).  So  let's  choose  a second  order  model: 

x(k)  = x(k-l)  + 02  x(k-2)  + u(k) 


where : 


1 p cos  9 

p COS  9 

p cos  9 1 

O' 2. 

2 

p cos  9 

Therefore : 

2 

pcos9(l-p  cos  26) 

“l  ■ 2 2T 

1 - p cos  9 

2 . 2„ 

a > -P  sin  9 

^ 1 - p^  cos^9 
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and  by  using  2-144:  The  variance  of  u(k,il)  is 


E{u^ (k, 
for : 


b) 

x(k) 

‘ for; 

I 

1 


*-)}  = Q = 1-a^^p  COS0  - U2P^  cos^9 

P = 0.96 

0 = 8“ 


= 1.127 


ct2  = -0.1854 


Q 


0.0911 


.variance  of 
Vodeling  error 


Eq.  2-71,  2-72,  lead  to  a model: 

2 

= p x(k-l)  - pcos  9 x(k-2)  + u (k) 

cl 

= aj^'x(k-l)  + a2*x(k-2)  + u^(k) 

p = 0.96 

9 = 8“ 
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0,9216 


a^'  = -0.95 


r.  T cn  ^ .Variance  of  , 
Q - 0.150  — ^ ^modeling  error 


= /Q’  = 0.387 


Q'  is  found  from:  Q.'  = (l-p^)(l  + p^).  Table  2 

summarizes  the  results. 


Table  2: 


R (n)  = p'  ' cos 

XX 

(9n) 

p = 0.96 

O) 

II 

00 

0 

Optimal  estimation  approach 

Filter  Response  method 

One  can  choose  to  represent 
x(k),  theoretically,  by  a 
set  of  infinite  previous 
points  (the  order  of  the 
difference  equation  can 
be  infinite) . 


X (k)  =a^^x  (k-1)  +U2X  (k-2 ) +u  (k) 


2 

_pcosa(l-p  cos2a)  _ i 
^ l-p‘^cos"'a 


Inherently  from  the  model- 
ing process  it  follows  that 
the  difference  equation  that 
expresses  the  process  is 
of  second-order. 


X (k)=aT 'x (k-1) +a_ 'x (k-2) +u_ (k; 

X ^ a 


I 


2 2 
-p^  sin^  9 


= -0.1854 


1 - p cos  ct 


var{u}  * modeling  error 

= v/l-Oj^pcose  - a2P^cos29 

- 0.302 


a^'  = p = 0.9216 


02'  = - P cos  9 =-0 . 95 


var{u}  = modeling  error 
= /(l-p^)  (1+p^) 


0.387 
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One  can  see  that  the  modeling  error  in  the  optimal  solution 
is  not  much  smaller  than  in  the  Filter  Response  Method 
(0.302  compared  with  0.387).  Therefore,  in  the  sense  of 
optimality,  the  model  in  Eq.  2-67,  2-68  is  "almost  good" 
as  the  optimal  solution.  It  seems  that  the  "nature"  of 
the  process  ~ P ^ ^ cos  (0n)  is  to  be  modeled  by 

Eq.  2-67,  2-68.  These  equations  are  simpler,  and  easy  to 
handle. 

Another  point  is  that  the  method  of  Section  F 
[Linear  Filter  Response]  leads  to  a representation  of  x{k,l) 
by  a finite  number  of  points.  Using  the  optimal  estimation 
approach,  one  can  find  a model  by  using  three,  four  and 
more  points,  previous  to  x(k,Z).  The  coefficients  of  the 
far-away  points  are  not  zero.  In  other  words  it  seems  like 
the  "nature"  of  this  case  is  not  to  be  presented  by  orthogonal 
projections. 

Finally,  these  models  are  used  in  recursive  estimation. 
In  this  application,  anyway,  the  whole  set  of  previous  data 
is  used  to  estimate  a point  (in  a recursive  way) . There- 
fore, both  models  lead  to  an  optimal-estimation-solution. 

We  can  summarize,  as  a consequence  of  this  example.  In 
the  sense  of  optimality,  the  Linear  Filter  Response  method 
is  worse  than  the  orthogonality  principle  method.  But  it 
has  other  obvious  advantages  that  make  it  a good  candidate 
for  signal  processing  and  especially  image  processing. 
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I . SUMMARY 

Note:  It  is  recommended  that  Example  5 in  Section  H 

be  studied  until  the  end  of  that  section, 
before  reading  Section  I. 

1.  This  chapter  has  developed  a state-variable  representation 
of  two  dimensional  random  fields. 

2.  The  concept  of  the  modeling  is  primarily  based  on  passing 
white  noise  through  a linear  filter.  The  basic  equations 
that  were  used  is  Eq.  2-40  and  2-41.  This  concept  was 
explained  in  Section  D,  and  then  was  used  in  Sections 

E,  F.  It  was  called  "filter  response  method".  The 
method  was  applied  to  three  types  of  Markov  fields: 

a)  The  first-order  type,  where: 

n / \ 2 Ini  I ml 

Rxx(n,m)  = ^ Pi  P2 

b)  The  "band  limited"  type,  where  the  spectrum  of  the 
signal  is  limited  in  a certain  band: 

Rxx(n,m)  ~ I " ^ cos  (9j^n)  ^ cos  (02ni) 

c)  The  second  order  type,  where: 

Rxx(n,m)  = (1+0^  ln|)P2'"'l  (1+02^1) 

= Unpj^  02  = J-np2 
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3. 


In  all  cases  it  was  shown  that  the  random  fields  can 


be  expressed  in  a form  that  was  recently  suggested  by 
Roesser  [Ref.  12]. 

M(k+l,Jl)  A2  M(k,Z) 

= + u (k,2,) 

N(k,£+1)  A^  N(k,Jl) 

The  nature  of  models  with  such  a structure  was  discussed 
in  Section  G. 

4.  The  "filter  response  method"  of  modeling  as  it  was 
discussed  in  Section  F has  two  disadvantages: 

a)  The  difference  equation  which  is  used  to  represent 
the  field  (i.e.;  the  linear  filter)  turns  out, 
always , to  be  a causal  equation.  But  in  image 
processing,  there  might  be  reasons  to  use  non- 
causal  filters.  It  is  obvious  that  representing 

a random  field  by  a non-causal  filter  will  be 
better  than  a causal  filter 

b)  This  method  is  limited  to  separable  autocorrelation 
functions  (smaller  modeling  error) . 

5.  The  modeling  method  by  the  orthogonal  principle  has 
its  disadvantages: 

a)  The  weighting  factors  have  no  simple  expressions. 

b)  The  number  of  neighbors  one  has  to  use,  is  theoretically 
infinite  for  some  types  of  autocorrelation  functions. 
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/ P 


c)  The  calculation  of  the  coefficients  required 


solving  many  algebraic  equations,  especially 
in  the  two-dimensional  case.  The  number  of 
equations  is  equal  to  the  number  of  coefficients, 
d)  Perhaps,  the  greatest  disadvantage  in  using  this 
method  for  modeling  two  dimensional  fields  is 
the  fact  that  it  is  difficult  to  represent  the 
two  dimensional  difference  equation  that 
describes  the  field  in  the  "state-space- 
structure"  [Ref.  12]. 

The  "Filter  Response"  method  does  not  suffer  from  these 
disadvantages . 

6.  Although  there  are  two  disadvantages  mentioned  above, 
it  seems  that  the  "Filter  Response  Method"  is  a good 
candidate  for  modeling  random  fields. 

7.  Table  3 summarizes  the  properties  of  the  two  methods: 
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Table  3:  Summary  of  Two  Methods  of  Modeling 


optiiral  estimation 
approach 


Possible 


"filter  response" 
approcich 


Ihe  model  is  not 
necessarily  optimal, 
and  the  modeling 
error  is  not  always 
white  noise 


Impossible 


see  Example 
5 - Section  H 


see  Example 
4 - Section  H 


Handling  of 
non-separable 
autocorrelation 
functions 


Possible 


Impossible 


see  Example 
3 - Section  H 


Gonplication 

of: 

- calculation 

- model  form 

Calculation  - com- 
plicated (many 
algebraic  equations) 

Model  form  - in  sate 
caises  complicated. 

Calculations  are 
sirtple  and  the 
model  has  a 
simple  structure 

see  Exanple 
5 - Section  H 

Deriving  a 
form  of  the 
model  by 
using  a 
state  vector 

Difficult  in  most 
cases  [even  for 
separable 
autocorrelation 
functions] 

Simple 

see  Example 
5 - Section  H 

Is  the 
recursive 
equation  that 
describes  the 
process  (model 
finite? 

Hieoretically,  for 
most  causes,  the 
recursive  equatiai 
is  infinite 

Yes 

see  Exanple 
5 - Section  H 

How  nany 
points  to 
include  in 
the  model 

'Iheoreticad.ly  - 
infinite. 

Practically  - enough 
points  to  make  the 
modeling  error  small 

The  order  of  the 
recursive  equation 
that  describes  the 
process  is  an 
inherent  result  of 
the  modeling  procedure 

see  Exanple 
5 - Section  H 
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III.  REVIEW  OF  ESTIMATION  THEORY  CONCEPTS 
A.  INTRODUCTION 

1.  Definition ; Stochastic  estimation  is  the  operation 

of  assigning  a value  to  an  unknown  system  state  or  parameters, 
based  upon  noise-corrupted  observation  involving  some  function 
of  the  state  or  parauneters. 

2.  For  example,  consider  a random  field  whose  model 
is  the  f irst-order-Markov: 

x(k+l,2,  + l)  = p^x(k+l,Jl)  -Pj^p2x(k,il)  + p2x(k,il+l)  + u(k,£) 

3-1 

where  u{k,£)  is  white  noise; 

{0  if  k?^ior£?^j 

3-2 

Q if  k=ior£=j 

The  observation  of  the  field  is  corrupted  by  noise:  y(k,£) 

is  the  observation: 

y(k,£)  = x(k,£)  + v(k,£).  3-3 

where : 


122 


0 if  m ^ i 


or  n / j 


E{v  (k,  J.)  V (i,  j ) } = / 3-4 

(_R  if  n = i and  n = j 

The  problem  is  to  find  the  "best"  estimate  for  x(k,l)  and 

A 

the  so-called  x(k,i),  from  a set  of  noisy  measurements 
y(i,j),  i = 1,  ...  k,  j = 1,  ...  2,. 

3.  In  Section  B,  the  criteria  for  an  estimator  to  be 
the  "best"  (optimal)  are  defined. 

Section  C defines  the  "Linear  Estimator". 

Section  D defines  the  most  common  estimator:  the 

linear  estimator  with  the  quadratic-cost  criteria.  The 
recursive  and  non-recursive  estimators  are  discussed. 

The  orthogonality  principle  is  also  explained. 

This  chapter  explains  the  concepts  of  optimal 
estimation.  Some  of  these  concepts  will  be  used  in  the 
next  chapter  to  solve  a few  particular  problems. 

B.  ESTIMATION  CONCEPTS  AND  CRITERIA 

In  this  section  "optimal  criteria"  will  be  defined. 

1.  Baysian  Estimation 

The  problem  is  to  estimate  a state  X,  from  a set 
of  measurements  Y. 

P(x)  - is  called  "a  priori"  probability  density 
function. 
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P(x/y)  - is  called  "a  posteriori"  probability 
density  function.  Knowing  the  density  function  P(x,y), 
various  types  of  estimators  can  be  determined.  Fig.  29 
will  be  used  to  define  them: 


Fig.  29:  The  Various  Types  of  Estimators 

a)  The  most  probable  estimate  is  (the  x 
that  is  most  likely  to  happen) . 

b)  The  estimator  is  X2»  the  solution  of 
minimizing  a conditional  mean.  The  problem  is  to  find 

00 

min  / 4)(x-x)  P(x/y)  dx  3- 


<)|  (x  - x) 

is 

called  the 

cost  function.  The  following  are 

typical 

e.Teunples . 

1) 

(x  - X)  = 

|x-  x|  = absolute  value. 

2) 

I)  (x  - x)  = 

(x  - x)  = square  error. 
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Both  types  are  the  so-called  equal  risk  estimators.  The 
probability  of  x being  larger  or  smaller  than  X2  is  the  same. 

The  quadratic  form  is  the  most  commonly  used  cost 
function. 

A 

c)  The  minimax  estimate  is  x^  ~ the  estimate  that 
minimizes  the  maximum  probability  of  error,  | x - x^ | . This 
is  simply  the  median. 

The  Bayes  Rule: 

Bayes  rule  is: 


P(x/y) 


P (x,y)  ^ P (y/x)P (x) 

P(y)  P(y) 


i 


i 

P(x) 


a priori  density  function  of  x 


P(y)  = density  function  of  the  measurements 

P(x,y)  = joint  density  function. 

P(y/x)  = conditional  density  function  of  the 
measurements  y. 

2.  Maximum  Likelihood  Estimation 

The  fundamental  idea  is  to  define  a "likelihood 
function"  and  to  maximize  this  function. 

The  estimation  problem  assumes  the  form  of  the 
probability  density  function  of  y,  P(y) , is  known: 


P(y)  = P(y,x) 
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where : 


1 


i 


i 


X is  unknown  parameter 
y is  a set  of  measurements. 

Therefore,  one  can  ask  what  is  the  best  estimate  of  x based 
on  these  measurements?  To  find  the  answer  it  has  to  be 
found  what  is  the  set  of  parameters  x,  that  will  cause  the 
set  of  measurements  y = (y^^  ...  y^)  most  likely  to  occur. 
The  following  procedure  of  computation  is  done: 

a)  Set  up  a likelihood  function 

L(y,x)  = P(y/x) 

b)  Find  the  parameter  x which  maximizes  this 
likelihood  function.  Hence  the  solution  must  satisfy  the 
condition : 


3L 

9x 


0 , 


< 0 


C.  LINEAR  ESTIMATION 

Optimal  estimation  theory  is  the  subject  of  finding  an 
estimator  that  minimizes  a given  "cost  function" , or  that 
is  "most  likely"  to  occur.  It  is  very  common  to  add  an 
additional  requirement  to  the  estimator;  The  requirement 
of  linearity.  The  linearity  requirement  will  be  introduced 
here  in  two  ways:  Non- Recursive  and  Recursive. 
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1 


1.  Non-Recursive  Estimation 

If  Y is  a given  set  of  measurements,  then  the 
estimator  will  be  a linear  combination  of  this  measurement: 


St  = Z a.y.  3-6 

1 1 

ley 

As  a particular  exaunple,  assume. an  image  that  is  scanned  top 
to  bottom,  left  to  right,  and  that  a causal  estimator  is 
required.  In  that  case: 


Fig.  30: 


Definition  of  Region 
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The  problem  is  to  find  the  coefficient  a(p,q)  that  will 
make  the  estimator  optimal  in  the  sense  that  was  described 
in  the  previous  section. 

2 . Recursive  Estimation 

For  easy  implementation  it  would  be  desirable  to 

A 

express  x(k,i!,)  in  terms  of  previous  estimators.  For 
example: 


x(k,a)  = I d.  (,  (p,q)x(p,q)  + G,  y(k,i) 

(p,q)£D(k,i) 


3-8 


The  region  D(k,i)  is  seen  in  Fig.  22.  Now  the  problem  is 
to  find  the  d's  and  G(k,2.)  for  each  point  in  the  image, 
satisfying  an  optimal  solution  (with  a given  cost  function) . 
It  is  easy  to  see  that  3-7  could  be  expressed  in  the  form 
of  3-8  (the  a's  in  3-7  could  be  expressed  in  terms  of  the 
d' s and  the  G' s) . 

There  is  no  doubt  that  the  non-recursive  estimator 
is  the  best  that  could  be  done,  simply  because  it  includes 
the  whole  set  of  given  data  in  the  "best"  way.  The  recur- 
sive estimator  also  uses  the  whole  set  of  data  to  estimate 
a point  x(k,£).  It's  only  done  in  a recursive  way.  There- 
fore the  best  recursive  estimator  that  can  be  found  is  the 
one  in  which  the  "cost"  is  equal  to  the  non-recursive 


counterpart. 


Of  course  one  can  ask  if  such  a recursive 
estimator  (whose  cost  is  equal  to  the  non-recursive 
counterpart)  exist  ? The  answer  to  this  question  empha- 
sizes the  difference  between  one-dimensional  estimation 
and  two  dimensional  estimation. 

In  the  one  dimensional  case,  if  the  model  of  the 
process  is  linear  and  the  cost  has  a quadratic  form,  both 
ways  of  estimation  have  the  sair.e  cost  (same  variance  of 
error) . In  that  case  the  recursive  estimator  is  the  well 
known  Kalman  filter . 

In  the  two  dimensional  case,  under  the  same  condi- 
tions, there  is  no  way  to  find  a recursive  estimator  that 
is  as  good  as  the  non-recursive  estimator.  It  will  be 
proven  in  the  next  chapter.  On  the  other  hand  it  will  be 
shown  that  one  can  find  a recursive  estimator  whose  variance 
of  error  is  almost  as  low  as  the  variance  of  the  non-recursive 
counterpart.  Therefore  the  recursive  estimator  will  be  a 
good  candidate  for  estimation  of  two  dimensional  random 
fields . 

The  next  section  will  define  the  most  commonly 
used  estimator. 

D.  LINEAR  ESTIMATION  COMBINED  WITH  QUADRATIC  COST  FUNCTION 
1 . Definition  Of  The  Estimator 

Let  the  cost  function  (|)  (x  - x)  in  3-5  be  a quadratic 

Ak 

function  (or  "quadratic  form")  of  (x  - x) . The  optimal 
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estimate  is  obtained  by  minimizing  the  cost; 


E{  [x(k,H)  - x(k,il)  ]^}  3-9 

Now,  this  cost  is  combined  with  another  restriction:  The 

estimator  should  be  linear. 

2 . The  Non  Recursive  Case 

The  linear  estimator  is  given  in  3-6 

X (k, 2.)  = z y (p»q) 

(p,q)€f2^  (k,2) 

The  error: 

^(k,2)  = x(k,Z)  - x(k,2)  = x(k,2)  - ' E a y(p,q) 

(p,q)ci^j^(k,2) 

and  the  variance  of  the  error: 

E{^>^(k,2,)}  = E{[x(k,Jl)  - E a _y(p,q)]^} 

(P,q)«J7j_(k,£) 

Now,  differentiating  with  respect  to  a , p = 1,  ...  k, 

p * q 

q = 1,  ...  I,  and  letting  the  derivative  equal  to  zero,  the 
following  set  of  orthogonal  conditions  must  be  satisfied 
for  optimal  solution: 

f 

E{(x(k,i)  - x(k,2)]  y(p,q)}  =*  0 
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or 


E{[x(k,Jl)  - Z A y(k,j)}  = 0 


{i,j)zQ^{k,D  3-10 


This  is  the  "orthogonal  principle"  in  optimal  estimation. 

The  meaning  of  3-10  is  that  the  estimation  error  x(k,£.)  - k(k,£) 
is  uncorrelated  to  the  set  of  data  that  is  used  to  estimate 
the  value  x(k,2.) 

To  summarize;  An  optimal  estimator  will  be  only  the  one  that 
satisfies  the  orthogonality  principle. 

3 . The  Recursive  Case 

It  is  clear  that  the  optimal  recursive  estimator  will 
be  the  one  that  satisfies  the  orthogonality  principle.  The 
problem  is  that  it  is  impossible  to  find  a two-dimensional 
recursive  estimator  that  does  it.  In  the  next  chapter  it 
will  be  shown  that  such  an  estimator  does  not  exist. 

Therefore,  the  approach  in  this  thesis  is  to  define 
a structure  for  the  filter,  say  as  in  3-8  [it  will  be  shown 
that  the  definition  of  the  structure  itself  leads  to  a 
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conclusion  that  the  estimator  is  not  optimal] , and  then 
to  find  the  parameters  of  the  filter  to  minimize  the  variance 
of  the  error. 
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IV.  RECURSIVE  ESTIMATION  OF 
TWO-DIMENSIONAL  FIELDS 


A.  INTRODUCTION 

1.  Problem  Definition 

In  this  chapter  the  general  description  of  the  problem 

will  be; 

a)  x(k,Z)  is  a stationary  random  field.  The 
autocorrelation  function  R (n,m)  is  known. 

b)  From  the  autocorrelation-function  one  can  find 
a dynamic  model  (state  variable  representation)  of  the 
field. 

c)  The  measurement  of  the  field  includes  noise: 

y (k,i,)  = X (k, «,)  + V (k,  i) 

x(k,il)  is  the  correlated  image.  v(k,ii)  is  white  noise. 
y(k,i)  is  the  measurement. 

d)  The  problem  is  to  estimate  x(k,il)  from  the 
measurement,  using  a linear  recursive  filter  that  minimizes 
the  quadratic  form 

E{  (x  - x) 

X is  the  estimate  of  x. 
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2.  General  Remarks  About  The  Solution 


1)  As  a consequence  of  Chpt.  Ill  one  thing  has  to 

be  kept  in  mind:  That  the  optimal  solution  is  defined  only 

in  a non-recursive  form,  by  Eq.  3-7,  to  minimize  the  mean 
of  the  square  error  (3-9) . This  definition  leads  to  the 
orthogonal  principle.  Therefore  the  optimal  recursive  filter 
will  be  only  the  one  that  is  proved  to  satisfy  the  orthogonality 
principle  [the  estimation  error  is  uncorrelated  to  the  data 
that  is  used  to  find  the  estimator] . In  the  one-dimensional 
case,  a recursive  filter  that  satisfies  the  orthogonal 
principle  exists  [this  is  the  well  known  Kalman  Filter] . 

But  in  the  two-dimensional  case,  if  one  defines  a structure 
for  a filter  that  is  equivalent  to  the  one  dimensional 
counterpart,  the  orthogonal  principle  cannot  be  satisfied. 

2)  The  conclusion  is  that  it  is  impossible  to 
find  an  optimal  recursive  filter , and  we  have  to  look  for 
an  sub-optimal  solution  [in  the  two-dimensional  problem] . 

3)  In  a sub-optimal  solution,  the  way  to  work  is 
to  define  a reasonable  structure  for  the  filter,  and  then 
to  find  the  gain  in  the  filter  structure  to  minimize  the 
mean  of  the  square  error. 

3 . Previous  Work  in  Two-Dimensional  Recursive  Estimation 

There  are  two  earlier  works  in  this  subject  [Ref. 

10  and  13].  The  appr 'ach  in  these  algorithms  was  to  assume 
that  the  orthogonality  condition  is  satisfied  for  all 
points  (i,j)  where  i £ k,  j ^ 1,  and  then  to  find  the  gain 
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for  point  k,  I , ^ induction.  These  algorithms  are 
incorrect  because  the  assumption  of  the  induction  is 
incorrect:  the  orthogonal  principle  cannot  be  satisfied 

for  point  (k,  il)  in  a recursive  filter. 

4 . The  Primary  Difficulty  of  the  Solution 

That  is  Suggested  in  this  Thesis 

The  algorithm  that  is  suggested  in  this  thesis  is 
not  "clear"  from  difficulties.  As  in  the  Kalman  Filter,  the 
suboptimal  gain  is  found  by  using  a recursive  set  of 
equations  [these  equations  calculate  the  variance  of  error 
and  the  gain] . But  the  problem  was  that  it  was  impossible 
to  find  a complete  "closed"  set  of  recursive  equations,  and 
one  approximation  had  to  be  done.  In  general,  it  is 
dangerous  to  make  an  approximation  in  a recursive  algorithm; 
a small  error  in  one  step  can  go  through  an  "integration 
process",  and  the  solution  might  "blow  up"  to  an  incorrect 
solution.  The  approximation  that  was  done  in  this  thesis  was 
tested  carefully  and  was  found  to  lead  to  excellent  results. 

B.  A RECURSIVE  FILTER  FOR  THE  PROCESS  R^^(n,m)  = 

1.  Introduction 

First,  let's  summarize  the  details  from  previous  chapters 
about  this  process  [using  the  definition  of  system  coordinates 
in  Fig.  26] . 

1)  The  autocorrelation  function  is 

R^^(n,m)  = P2  (Eq.  2-78) 
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2)  The  "model"  of  this  process: 


x(k+l,<!-  + l)  = p^x(k+l,)l)  - p^P2x(k,i)  + p^xCkfll  + l)  + u{k,)l) 

= varfu(k,i,)}  = /q  = a ^ (1-p  ^ (1-p 

S 12 

(eq's.  2-81,  2-82).  x(k+l,2.+l)  depends  on  his  neighbors 

in  region  D(k+l,Jl+l).  [Fig.  27.] 

3)  Each  line  and  column  of  the  field  has  a one- 
dimensional model: 


variance  of 
modeling  error 


x(k+l)  = p2x(k)  + U2(k) 
x(£  + l)  = p^x(2.)  + 


for  any  I 
for  any  k 


var{u2}  = |/l  - P2^  = ^ 

var{uj^}  = l/l  - p^^  = 


(eqs.  2-51,  2-53,  2-54) 


2 . Statement  of  the  Problem 

Given:  a discrete,  random,  two  dimensional  process: 


x(k+l,H+l)  = pj^x(k+l,£)  - p^P2X(k,£)  + p2X(k,£+l)  + u(k,£) 


4-1 
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and  the  noisy  observation,  starting  at  = (1,1) 


y(k+l,2.+l)  = cx(k+l,?,+  i)  + v(k+l,Z+l) 


4-2 


where : 


E{u(k,  £)}=  E{v(k,£)}  = 0 for  all  (k,£) 


4-3a 


E(u(i,  j)  v(k,Jl)  } = 0 for  all  i,j,k,il 

j"  Q if  k = i and  *-  = j 

1^0  if  k ^ i or  Ji  j 


E(u(k,  2.)u(i,  j)  } = 


4-3b 


4-4 


E{x(k,Jl)v(i,  j)  } = 0 for  all  k,2,,i,j 


4-5 


E{x(k,0)  = E{x(0,£)}  ^ X for  all  k,l 


4-6a 


E{  Cx(k,0)  - X]  = E{[x(0,£)  -x]^}  = o ^ for  all  4-6b 

® k,2. 


r 


E{v(k,Jl)v(i,  j)  } 


R if  k = i and  *■  = j 


0 ifk/i  or  Z ^ 3 


It  is  convenient  to  define: 


4-7 


(^1  ”^1^2  ^2^ 


4-3 
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x(k+l,£) 


x(k,?,)  = 


x(k,  M 
x(k,«.+l) 


using  4-1,  4-9,  4-10: 


4-9 


1 

: x(k+l,£,  + l)  = p x(k,Jl)  + u(k,£.) 


4-10 


Problem:  Find  an  estimate  of  x(k,l),  denoted  x(k,^)  which 

is  a linear  function  of  all  observations  y(i,j), 
minimizing  the  quadratic-form 


E{x(k,£)  - x(k,Jl)} 


4-11 


The  estimate  should  be  done  in  a recursive  filter. 

3 . The  One-Step  Predictor 

In  the  one-dimensional  filter,  the  one-step  predictor 
is  defined: 


x(klk-l)  = x(k-l|k-l) 

x(k|k-l)  means  the  estimate  of  x(k)  given  k-1  points 
x(k-l|k-l)  means  the  estimate  of  x(k-l)  given  k-1  points. 

Therefore,  it  will  be  reasonable  to  define  the  one 
step  predictor  in  the  two-dimensional  case: 
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I X (k+1, i+1)  = p^k(k+l,i)  - p^P2x(k,£)  + p2x(k,Z-l) 

^ T-12 


where: 


X (k+1, i+1)  is  the  estimate  of  x(k+l,H+l)  by  using  the 
data  of  n^(k+l,2+l) 

k(k+l,C)  is  the  estimate  of  x(k+l,H)  by  using  data 

of 


ji(k,Jl)  is  the  estimate  of  x{k,l)  by  using  data  of 

^2 

k(k,il+l)  is  the  estimate  of  x(k,il+l)  by  using  data  of 

f22(k,i+l) 


Defining : 


^(k,il) 


k(k+l,l) 
X (k, «.) 
k(k,i+I) 


4-13 


and  by  using  4-9,  4-12,  4-13: 


3^^^^  (k+l,^+l)  » p k(k,il) 


4-14 


But  there  is  a fundamental  difference  between  the  one 
dimensional  and  the  two  dimensional  one-step  predictor.  This 
difference  is  emphasized  in  the  next  theorem. 
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Theorem;  1)  The  one-dimensional  one-step-predictor  is 
optimal  if  x(k)  is  optimal. 

2)  The  two  dimensional  "one-step  predictor"  is 

not  optimal,  even  if  k(k+l,il),  k(k,)i),  k(k,Jl+l) 
in  equation  4-12  are  optimal. 

Proof ; See  Appendix  A. 

Conclusion:  The  recursive  filter  cannot  be  optimal. 


4 . Estimator 

When  one  has  the  measurements  of  the  point  (k+l,il+l) 
then  the  estimator  k(k+l,J,+l)  will  be; 


k(k+l,^+l)  = (k+l,)l+l)  + G(k+l,H+l)  [y(k+l,£+D  - Cx^^^  (k+l,£+l)  ] 


4-15 


The  effect  of  the  GAIN,  G(k,£)yis  to  weight  the  correction 
term  [y(k,S,)  - Cx  (k+1,  £+1)  ] . 

Remark;  Because  the  Estimator  and  the  one  step  predictor  are 
recursive  in  their  nature,  it  is  easy  to  see  that  the  value 
x(k+l,£+l)  is  actually  formed  by  the  measurements  in 
^2  *-+i)  and  k (k+1 , £+1)  is  formed  by  the  measurement 

in  (k+1,  £+1)  . 
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(1) 

(k+l,Jl+l) 


(i,  j)enj^(k+l,z+l) 


X(k+l,a+l)  = 


^ 6.  .y(i, j) 

(i,  j)sf^2  (k+l,i+l) 


Let  us  define: 


F(k+1,£+1)  = 1 - C-G(k+l,i+l) 


4-16 


Using  4-16,  4-15: 


x(k+l,i+l)  = F(k+l,i+l)x'  ' (k+l,i+l)  +G(k+l,l+l)y(k+l,l+l) 


4-17 


Now,  because  the  one  step  predictor  cannot  be  optimal,  the 
estimator  is  also  not  optimal  (there  does  not  exist  a value 
G(k+l,J,+l)  to  make  x(k+l,il+l)  optimal).  Of  course  one  can 
minimize  the  variance  of  error  of  the  filter  that  is  defined 
in  4-12,  4-15  and  get  a suboptimal  estimator. 

5.  Estimation  Error 
Define : 


jEstimationj  4 6(k+l,i+l) 

Error 


I 


= x(k+l,^  + l)  - x(k+l,2,  + l)  4-18 

and : 

f (k+l,£) 

d(k,il)  = r(k,^)  4-19 

£{k,)l+l) 

substituting  4-1,  4-2,  4-9,  4-15,  4-19  into  4-18: 

<^(k+l,£+l)  = F(k+1,£+1)  p§(k,£)  - F(k+l,£+l)u(k,£) 

+ G(k+l,I'+l)v(k+l,£+l)  ) 4-20 

In  the  situation  at  hand,  the  estimate  is  a random  process 
and  so  is  the  estimation  error.  Next,  we  shall  discuss  the 
statistics  of  the  estimation  error. 

6.  Unbiased  Estimator 

Using  4-20,  let  us  compute  the  mean  of  »(k+l,£+l): 
E{^(k+1,£+1) } = F(k+l,£+l)E{p^(k,£) } - F (k+1, £+1) E{u (k, £) } 

+ G(k+l,£+l)E{v(k+l,£+l) } 

The  last  two  terms  in  the  above  equation  are  zero  (using 
4-3a) . Therefore: 
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E{^(k+1,«.+1)  } 


F(k+l,^l+l)E{p£(k,^)  } 


= F(k+1,£+1)  [p^6:(k+l,^)  -p^P2^(k,£)  +P2^(k,£-1)] 

Now:  we  noticed  that  the  estimates  of  points  x(0,£),  x(k,0) 

are  not  based  upon  measurements  (recall  that  measurements 
are  given  for  k,£  ^ 1).  These  are  arbitrary  choices. 

Therefore  let  us  choose: 


X (o, £) 


X (k, 0)  = X 


.mean  of  the  , 
^random  field ^ 


In  this  case: 


^(k,0)  = Jr  - x(k,0) 


€(0,h  = X - x(0,£) 


and  using  4-6a 

S(k,0)  = f(0,k)  = 0 

Because  4-20  is  a recursive  equation,  it  is  easy  to  see: 
E{5(k,£)}  = 0 for  all  k,£. 

This  is  the  case  of  an  unbiased  estimator. 
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7.  Variance  of  Estimation  Error 


Definitions : 


4-23 

4-24 

4-25 

4-26 


One  has  to  be  careful:  4-23  is  a definition  of  a scalar 

P (k+1 , J?,+l)  . 4-25  is  a definition  of  a matrix  P(k+l,Jl+l). 

The  sign  (k,  £)  is  the  transpose  of  (•)(]<, Z).  P(k+l,Z+l) 
is  the  variance  of  the  estimation  error  at  point  (k+l,£+l). 
P^^^  (k+l,£+l)  is  the  variance  of  the  "one  step  predictor" 
error  at  point  (k+l,£+l).  P(k,£)  is  a 3x3  symmetric  matrix. 

Three  of  its  values  are  the  variances  of  the  estimation 
errors  of  region  D(k+1,£+1).  Three  other  values  of  the 
matrix  are  correlations  between  the  estimation  erros  of 
points  in  region  D(k+1,£+1).  For  better  understanding 
see  also  Eq.  4-27.  P.  ,(k,£)  is  the  correlation  between 
the  estimation  errors  of  points  (k,£)  and  (k+i/£+j). 
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Remark : The  definition  in  Eq.  4-23  is  a special  case  of 


Eq.  4-26; 

P(k,Z)  = PQ^g(k,l). 

8.  Variance  Calculations 

Using  the  definitions  in  4-23,  4-24,  4-25,  4-26, 
one  can  derive  recursive  equations  as  follows: 

P(k,l)  calculation 


P(k,il)  = E{f(k,^)  ^"^(k,)!)} 


= E 


( &(k+l,il)6(k+l,a) 
8(k,i)8(k+l,^) 

, 6(k,)l+l)€(k+l,jl) 


) 


6(k+l,^)^:k,£) 

g(k,Jl)€(k,Jl) 

€(k,JH-l)^(k,Z) 


^(k+l,l)(f(k,J,+l) 

€(k,^)€(k,Jl+l) 

^(k,£+l)^(k,il+l) 


using  4-26,  4-23: 
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P (k+1 , il+1)  calculation 

P^^Nk+l,l+l)  = E{  [x^^^k+l,il+l)  - x(k+l,Jl+l)]^} 

= E{[px(k,i)  -px(k,^)  -u(k,i)]^} 

Because  the  term  in  the  rectangular  brackets  is  a scalar 
we  can  rewrite: 

P^^Nk+l,£-t-l)  = E{[p^(k,^)  - u(k,Jl)  ] [p;^  (k,Z)  -u(k,^)]'^} 

= E{p8>(k,il)6^(k,)l)p'^  - u(k,^)6'^p'^  - p^(k,i)u(k,^)  +u^(k,^)} 

The  expectation  of  the  two  middle  terms  in  the  last  expression 
is  zero  [see  Appendix  B] . Therefore, 

P^^Nk+l,£+l)  = pE^!!(k,^)f^(k,^)  }p'^  + E(u^  (k,^)  } 

By  using  4-4,  4-25: 

4-28 

P(k+l,£-*-l)  calculation 

By  substituting  4-20  in  4-23  and  noting  that 
^*^(k+l,£+l)  = ^[k+l,£+l),  since  it  is  a scalar: 
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P(k+1,£+1)  = E{e^(k+1,£+1) } 

= E{[F(k+l,£+l)p(f(k,£)  -F(k+l,£+l)u(k,£)  +G(k+l,£+l)v(k+l,£+l)  ] 

• [F(k+l,£+l)p^(k,£)  -F(k+l,£+l)u(k^  +G{k+l,£+l)v(k+l,£+l)  f} 

The  last  equation  includes  the  three  next  terms  that  are 
zero : 

E{  i^(k,  £ ) V (k+1,  £+1)  } = 0 (see  Appendix  B) 

E{$  (k,£)u  (k,£) } = 0 (see  Appendix  B) 

E{u (k, £) V (k+1, £+1) } = 0 (by  using  3-3b) 

Therefore : 

P(k+1,£+1)  = F^  (k+1,  £+l)p  ^(k,£  )€'^{k,  £)  p + F^  (k+l,£  + l)u^  (k,£) 

+ G^(k+l,£+l)v^(k+l.£+l) 

and  by  using  4-4,  4-25,  4-7,  4-28: 

P(k+1,£+1)  = F^(k+1,£+1)  P^^Nk+l,£+l)  + G^  (k+1^  £+!)■  R 

4-29 
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Pq  Q(k,Jl)  calculation 

From  4-26: 

^i+1,  j+l^’^'*'^  = E{(§(k,l)^(k+i+l,Z+j+l) 

substituting  <f(k+i,2.+j)  from  4-20; 

Pi^l^  j^.1  (k,il)  = E{^(k,i)F(k+i+l,Z+j  + l)p§(k+i,£+j) 

+ E{  ^(k,Jl)G(k+i+l,a+j+l)v(k+i+l,i+j  + l) 

- E{€.(k,Z)F(k-t-i+l,£4-j+l)u  (k+i,£+j)  } 

The  last  two  terms  are  zero  [see  Appendix  B] . Now,  by  using 
4-26,  4-20,  4-19  and  4-9,  the  last  term  becomes: 


’i+1,  j + l(k,£)  = F(k+i+l,£+j+l)  (k,i)  - p^p^P^^  j (k,Jl) 


Substituting  in  the  last  equation  the  obvious  identities: 


P (k,Jl)  = P (k-i,^-j) 

/ “J  1 f J 

P (k,Jl)  = P (k-i,^+j) 
» j 1 / ” j 
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= F(k+1,^)  -p^P2Pg^j^(k,l+l) +P2Pg^Q(k,i)] 


calculation  of  P ^ , (k,?,) 

It  is  not  possible  to  find  a recursive  equation 
of  P_j^  ^{k,Jl)  that  includes  just  previous  values  of 
P(i,j),  Pq  ^(i,j),  P^  Q(i,j).  Trying  to  evaluate  P_^  j^(k,Jl) 
leads  to  expressing  P , , (k,Jl)  in  terms  of  correlation 
of  points  farther  removed  from  the  region  D(k,i).  This 
problem  is  solved  by  an  approximation . But  first  let's 
see  some  special  cases: 

a. 


= [x(0,D  - x(0,l)  ] tx(l,0)  - x(0,l)] 


and  using  2-78 
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b.  _j_(0,i) 

p,  . (0,i)  = E{^(1,  -l)^’(0,^) 

= E{%(l,Jl-l)  [x(0,«.)  -x(0,Z)]}  4-33 

Using  2-53,  2-54,  and  4-21: 

Pj_  _l(0,a)  = E{^(1,)1-1)  [p^x(0, «.-!)+ (0,^-1)  - x] 

Now,  using  the  same  procedure  as  in  Appendix  B it  can  be 
shown : 

E{^(l,^-l)u^(0,^-l) } = 0 

and  using  Eq.  4-22: 

E{8(l,il-1)  - x}  = 0 

Therefore: 

_^(0,«.)  = E{^(l,^-l)  pj_x  (0, 2.-1)  } 

= E{^(l,^-l)6,(0,i-l)'P^}  4-34 
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using  4-26: 


P(W,1-1) 


4-35 


c.  _j^(l,k) 

In  the  same  procedure: 


P2Po,i(k,0) 


P (k,61 

4-36 


d.  In  The  General  Case 

Comparing  Eq.  4-1  and  4-20  one  can  see  that  the 
error  is  a random  field  that  has  the  same  nature  as  the 
original  field  x(k,il).  Therefore  one  can  assume  that  the 
autocorrelation  function  of  the  error  has  a similar  form  to 

4-1. 

Rg^(n,m)  = 4-37 

This  autocorrelation  function  is  correct  only  in  the 
steady-state , when  F(k,2.)  is  constant.  Therefore  Eq. 

4-37  is  an  approximation  that  is  correct  in  steady  state. 

From  4-37  it  is  seen  that  in  steady  state: 

^1,0^^'^^  ’ P(k,i).p2  4-38a 


4-38b 


From  4028a, b,c: 

4-39 

Note  that  this  equation  is  accurate  for  first  line  and 
colume  (e.g.  4-35,  4-36). 

9 . Minimum  Variance  Condition 

We  have  found  the  condition  for  an  unbiased  estimator. 
This  condition  doesn't  require  any  constraint  upon  the  gain 
G(k,?).  Thus  we  wish  to  minimize  the  value  of  P(k+1,/+1), 
with  respect  to  G(k+1,£+1).  Differentiating  the  expression 
for  P(k+l,)l+l),  set  the  result  to  zero  and  solve  for 
G(k+1,?,+1).  From  4-29: 


° " 3G(k,£) 
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► 


And  substituting  4-40  into  4-29: 


P(k,Jl)  = (1  -G{k,l)  (k,^) 


4-41 


10.  Initial  Condition 

To  start  the  recursive  process  of  Gain  calculation, 
we  need  the  next  initial  conditions: 

A 

1.  x(0,)l)  0 < < N 

2.  x(k,«,)  0 < k < N 


and  then  we  are  able  to  calculate: 


3.  P(0,JI) 

4.  P(k,0) 

5.  Po^i(0,^) 

6.  Pj^^g(k,0) 

7. 


0 £ Jl  ^ N 

0 ^ k ^ N 

0 ± I ± N-1 

0 < k < N-1 


Fig.  30  shows  what  is  meant  by  'initial  condition." 
For  an  unbiased  estimator: 


x(0,)l)  » x(k,0)  = X 
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By  using  4-16,  4-21,  4-23; 

P(0,ii)  = E{£^(0,«,)}  = E{  (x(0,Z)  - x(0,)l)  )^} 

= E{  [x  - x(0,£) 

by  using  4-6b; 

P(0,i)  = oj 

- ■■■  S 

and  similarly  for  P(k,  0): 

P(k,0)  = 4-43 

For  • 

Po,i(0,£)  = E{e(0,£)e(0,£+1) 

ys  /V 

= E{  (x(0,£)  - x(0,£)  ) (x(0,£+l)  - x(0,£+l)  ) 
= E{  (X  - x(0,  £)  ) (X  - x(0,£)  ) 

PQ,l(0,t)  - 0;  cj  4-44 

and  in  the  s£une  form: 
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4-45 


Pl,.l(0,l)  = p^P2 


4-46 


11.  Summary 

The  proposed  filter  for  the  process  that  is  defined 
in  4-1,  4-2  is  given  in  4-12  and  4-15.  The  recursive 
equation  for  GAIN  calculations  are 


G(k,Z)  = P (k,  ^)C[C^P  i)  + R]“^ 


= p P(k-l,Jl-l)  p'^  + Q 


P(k,il)  =.  (1  - G(k,i)  ) P^^^  (k,^) 


where : 


P(k+l,i) 

P(k,Jl)  = 1 P^^g(k,l) 


P(k,Jl) 


Pl^-l(k,£+l)  Pg^^(k,£) 


P _ (k,£+l); 

Po,i(k,^)  ; 

P(k,£+1)  / 


(Pi  -P1P2 


'2' 


and: 


155 


I 


= F(k,i)  - p^P2Po,l^^~^'^~^^  +P2P(k-l,^) 


Pg^j^(k,£-1)  = F(k,)l)  [pj^P(k,Jl-l)  - p^P2Pj^^Q(k-l,i-l)  + 


P.  .(k,£)  P-  , (k,l) 
= ' P(k,t)°' 


Actually  it  is  impossible  to  find  a complete  set  of  recursive 
equations.  The  equation  for  P^^  _j^(k,Jl)  is  only  an  approxima- 
tion for  the  case  where  1^2  and  k ^ 1.  Therefore  we  have 
an  algorithm  that  is  correct  for  the  first  line  and  column 
and  converges  to  a correct  suboptimal  solution  in  steady 
state. 

Fig.  35  helps  to  understand  the  procedure  of  the 
recursive  calculation.  Assume,  we  have  finished  the  calcu- 
lations ,for  point  (k+l,H)  and  we  wish  to  calculate  the  gain 
G(k+1,£+1).  In  order  to  do  this  we  need  to  know  P (k+1 , £+1) 
[variance  of  one  step  predictor] . Now: 

The  value  P^^^  (k+l,i!.+l)  depends  on  six  values 
that  are  seen  on  the  boundaries  of  triangle  I in  Fig.  29; 


P(k,£+1),  P(k,£),  P(k+1,£),  Pq  J^(k,^),  P^^  Q(k,H),  (k,  ^+1) 


- The  value  P-  , (k,i)  depends  on  three  values  that 
u , 1 

create  the  boundaries  of  triangle  II  in  Fig.  29; 
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Po^l(k,H-l)  , P(k,Jl)  , 

- The  value  qOc^)  depends  on  the  three  values 
that  create  the  boundaries  of  triangle  III  in  Fig.  29: 

Pl^_l(k,£),  PQ^j^(k,Jl-l)  , PQ^Q(k,^). 

The  value  of  P^^  _j^(k,£+l)  is  given  in  Eq.  4-39.  As  one 
moves  along  a line,  this  structure  is  also  moving. 

11.  Results 

Figures  36,  37,  38  show  typical  results.  Steady 
state  is  reached  after  4-5  points.  The  steady  state  for 
the  first  line  and  column  is  higher  than  for  points  in  the 
middle  of  the  field.  The  reason  for  this  fact  is  that 
points  in  the  first  line  are  estimated  only  by  their  neighbors 
to  the  left.  It  was  found  that  the  results  for  the  first 
line  and  column  are  exactly  the  same  as  one  finds  from  the 
one-dimensional  Kalman  Filter,  using  the  model  in  Eq.  2-53, 
2-54. 

12 . Checking  Of  The  Results 

Two  things  must  be  checked:  First,  it  should  be 

checked  if  the  algorithm  developed  above  is  correct.  Recall 
that  it  was  impossible  to  develop  a complete  accurate 
algorithm  and  an  approximation  had  to  be  used.  In  any 
recursive  calculation  an  approximation  can  "blow  up"  the 
results.  Therefore  it  is  very  important  to  check  if  this  is 
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Fig.  31:  Definitio 

System 
Coordinat 
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K 
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Fig.  33:  Definition  of  Regions  (k+1 , 2,+l)  , ^2 


Fig.  34:  Initial  Conditions.  A black  point  in  a cell 

(i,j)  represents  the  initial  value  x(i,j)  and 
the  value  P(i,j).  A line  between  two  black 
points  represents  the  correlation  between  the 

estimation  error  of  the  two  pixels:  P.  .(k,)!,). 

!■ » 3 


GAIN  OF  FILTER 


not  the  case  here.  It  was  found  that  the  results  for  the 


first  line  and  column  are  the  same  as  the  one  dimensional 
Kalman  filter.  Now,  we  want  to  check  the  steady  state, 
after  MxM  points.  In  order  to  solve  this  problem,  the 
next  procedure  was  done: 

1)  An  image,  of  size  128  x 128  with  autocorrelation 

function 


R (n,m) 

XX 


m 


was  created. 

2)  White  noise  with  variance  /R  was  added  to 
the  correlated  image. 

3)  The  filter  developed  above  was  used  in  order 

to  estimate  the  correlated  process.  The  gain  of  the  filter 
was  found  by  substituting  d^,  R^into  the  recursive 

algorithm  of  gain  calculation. 

4)  The  variance  of  error  of  the  simulation  was 
compared  with  the  theoretical  variance  of  error. 

This  procedure  is  described  in  Fig.  39. 

The  basic  assumption  here  is  that  if  the  simulated 
variance  of  error  is  similar  to  the  theoretical  variance  of 
error,  then  the  algorithm  is  correct.  A mathematical 
justification  for  this  assumption  was  not  proved.  Still 
one  can  ask  if  the  mean  square  error  is  minimum.  But  the 
theoretical  variance  of  error  is  an  inherent  result  of 
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gain  calculation,  and  part  of  the  algorithm.  Therefore 
if  it  is  equal  to  the  simulated  result,  it  is  reasonable 
to  assume  that  the  algorithm  is  correct. 


’‘(k.l) 


SIMULATION  OF  AN  IMAGE 


Fig.  39:  Experiment  to  Check  the  Algorithm  for  Estimation 

by  Using  a Simulated  Image 

The  second  question  arises  from  the  knowledge  that 
the  two-dimensional  filter  is  not  optimal.  Therefore  it  is 
imp'  rtant  to  compare  the  optimal  non-recursive  estimation 
error  to  the  sub-optimal  recursive  filter. 

Results : 

Fig.  40  shows  three  types  of  results: 

optimal  non-r  cursive  estimation  error  (in 
steady  state) . 

theoretical  steady  state  variance  of  error  of 
the  recursive  filter. 
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experimental  steady  state  variance  of  error  of 
the  recursive  filter. 

2 

This  experiment  was  carried  out  for  = p2  = 0 . 96 , o =1. 
The  variance  of  random  noise  varies  from  0.1  to  0.8. 

It  was  found  excellent  coincidence  (less  than 
5%)  between  the  theoretical  estimation  error  and  the 
simulated  error.  The  recursive  filter  shows  variance  of 
error  not  bigger  than  2%  than  the  optimal  non-recursive. 
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theoretical  recursive  filter 

experimental  results  using  the  recursive  filter 
optimal,  non-recursive  estimation  (theoretical) 


V.  APPLICATIONS:  DETECTION  OF  TARGETS 

IN  PRESENCE  OF  CORRELATED  NOISE 

A.  INTRODUCTION 

Assume  a problem  of  detecting  an  infrared  target.  In 
this  case  one  has  interest  in  the  location  of  the  target 
and  the  intensity  of  the  target.  In  infrared  detection 
the  intensity  has  a special  importance,  because  the  target 
intensity  is  proportional  to  the  temperature,  and  by 
knowing  the  intensity  one  can  conclude  whether  the  target 
is  a missile,  or  an  aircraft,  etc.  The  target  detection  can 
be  interfered  by  two  types  of  noise: 

1.  Correlated  Noise.  An  example  for. correlated  noise 
is  clouds. 

2.  White  Noise.  An  exaunple  for  white  noise  is 
measurement  noise. 

Now,  the  problem  is  to  find  the  location  and  intensity  of 
that  target. 

B.  STATEMENT  OF  THE  PROBLEM 

In  this  chapter  the  next  problem  is  solved: 

Given:  A set  of  measurements  of  a two  dimensional  field, 

y(k,l).  These  measurements  are  composed  of  three 
types  of  signals: 

1)  A target  that  has  an  intensity  function  T(k,il). 

2)  The  target  is  corrupted  in  a correlated  background 
x(k,H).  The  autocorrelation  function  is,  for 
exeunple : 
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3)  White  noise,  vCk.,1),  that  is  uncorrelated  to 
the  background  or  target,  so  that: 

y(k,il)  = x(k,i)  + T(k,«,)  + v(k,£)  5-2 

Find : 

1)  The  location  of  the  target. 

2)  The  intensity  of  the  target. 

C.  BASIC  CONCEPT 

Assuming  the  autocorrelation  of  the  background  is  known, 
one  can  build  an  estimator  for  the  background  [see  Chpt.  IV] . 
The  filter  is  applied  with  steady  state  gain.  We  use  the 
"one  step  predictor"  and  compare  it  with  the  actual 
measurement. 

After  subtracting  the  estimated  background  from  the 
measured  image,  the  residual  image  includes  mainly  the 
noise  and  the  target.  If  rhe  target's  intensity  is  higher 
enough  from  the  noise,  the  detection  can  be  dona  in  the 
residual  image. 

A further  improvement  can  be  achieved  if  one  knows  the 
shape  of  the  target.  In  that  case  classical  methods  like 
matched  filters  will  increase  the  processing  gain.  Figure 
41  shows  the  basic  concept. 
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Fig.  41:  Basic  Concept  of  Target  Detection 

D.  DETECTION  OF  A POINT  TARGET 
1.  The  Two  Hypotheses 

From  Fig.  41  it  is  seen  that  the  "residual 
image"  is  composed  of  white  noise  and  the  target.  As  a 
consequence  one  can  see  that  a point  target  is  detectable 
with  a low  probability-of-error  only  if  its  intensity  is 
significantly  higher  from  the  white  noise,  say  three  times 
more  than  the  variance  of  the  noise. 

The  point  target  detection  is  done  as  follows: 
When  arriving  at  a new  point,  say  (k+l,Z+l),  during  the 
scanning  process,  the  situation  at  hand  is  that  there  are 
two  estimators  of  the  seune  point: 

^11) 

One  - is  the  "one  step  predictor",  x (k,£), 
that  does  not  include  the  measurement  at  the  point. 
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The  Second  is  the  measurement  itself,  y(k+l,2,+  l). 
Now,  if  y(k+l,!l+l)  is  significantly  higher  than  the  "one 
step  predictor",  it  is  most  likely  due  to  a target.  The 
target  detection  is,  therefore,  a result  of  comparing  the 
measurement  at  a point  to  the  prediction  of  the  value  at 
that  point,  by  using  all  of  the  neighbors  near  the  point. 

It  is  clear  that  this  difference  can  be  due  to  the  measurement 
noise  v(k,il).  Therefore  a reasonably  good  detection  will 
be  if  the  target  is  three  times  higher  than  the  variance 
of  the  noise. 

Now,  the  procedure  will  be  explained  mathematically: 
Assume  a specific  problem: 

1)  The  autocorrelation  function  of  the 

background : 


2)  The  measurement  noise  has  a Gaussian,  white 
noise  statistic: 


Rvv(n^"') 


l^R  if  n = m=  0 

[ 0 if  n ^ 0 or  m ^ 0 


3)  The  target  has  a Gaussian  probability 
density  function  with  mean  T^  and  variance  /R^  . 
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Note : 


It  is  more  reasonable  to  assume  a Rayleigh  distribution. 
The  Gaussian  assumption  is  done  because  it  simplifies 
the  discussion. 


Fig.  42;  Probability  Density  Function 
Of  the  Target 


The  estimation  of  the  background  is  done  by 
using  the  filter  that  is  described  in  Chpt.  Iv.  Also  assume 
that  the  procedure  of  estimation/detection  is  finished  for 
point  (k+l,H)  and  now  it  is  repeated  for  point  (k+l,Jl+l). 

Recall : 

A / 1 ^ • 

X (k+l,il+l)  = one  step  predictor. 

x(k+l,Z+l)  = estimator  of  point  fk+1,  ^+V. 
using  Eq.  4-15,  assuming  C = 1; 


i 


! 

I 

I 

i 
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3.  The  variance  of  the  one  step  predictor 


^ p(l) 

4.  ~ represents  hypothesis  no.  1 

presence  of  the  target  at  point  (k+l,2.+l). 

^2  - represents  hypothesis  no.  2 
absence  of  the  target  at  point  (k+1,2,+1). 

Due  to  these  hypotheses : 

r 

y(k+l,2,+l)  = x(k+l,<l+l)  + v(k+l,)l  + l)  + < 


and : 


r(k+l,£+l) 


= -y(k+l,£+l) 


+ x^'^(k+l,  £+1) 


= ,fe*^\k+l,£+l) 


+ V (k+1, i+1)  + 


L° 


5-4 

, the 
, the 


T when 

0 when  H2 
5-5 


when 

when  H2 
5-6 
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The  Processing  Gain 
Definition ; 

,processinq^  _ c i 

^ gain  ' (S/N) 

In  our  case:  The  intensity  of  the  target,  T,  as 

the  residual  image  is  equal  to  the  intensity  at  the  input 
(at  the  measurement  device) . 


out 


S . 
in 


and  therefore: 


-processing. 

gain 


N ^ 
out 

N. 

in 


The  variance  of  the  input  noise: 


N.  = /C  + R 
in 


a = variance  of  the  correlated  field. 
/R  = variance  of  the  white  noise. 


The  variance  of  the  output  noise: 


N ^ + R 

out 
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variance  of  "one  step  predictor"  error 


(1)  = 
P 


5-8 


The  Threshold  Device 

Finally,  one  has  to  define  two  regions  on  the 
r(k+l,2,  + l)  axis: 

- a region  of  values  where  the  decision  will  be  for 
presence  of  a target. 

- a region  of  values  where  the  decision  will  be 
for  absence  of  a target. 

The  input  to  the  threshold  device  is  the  "one  step  predictor" 

Case  1:  The  mean  T is  not  known. 

o 

-In  this  case  the  threshold  has  to  refer  only  to  the  variance 
of  r when  H2  occurs.  A reasonable  threshold  is,  for 
excunple 


threshold 


3 
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Case  2:  The  mean  T is  known. 

— o 

In  this  case  the  threshold  has  to  refer  to  T . In  this 

o 

case  one  can  use  a "window  threshold".  See  Fig.  45. 


probability  dansity  function 


WINDOW 


THRESHOLD 


▼ 

PROBABILITY 
OF  "miss" 


Fig.  45:  Threshold  for  Case  2 (T^  is  Known). 


E.  DETECTION  OF  LINES 
1.  The  Problem 

There  are  many  algorithms  to  detect  lines.  Most  of 
them  suffer  from  one  disadvantage:  that  to  detect  a line, 

an  a-priori  assumption  must  be  about  the  direction  of  the 
line.  Therefore,  in  order  to  detect  lines  in  several 
directions,  the  image  has  to  be  scanned  several  times  or, 
during  one  scan  many  calculations  must  be  done  for  each  point 


6 


[f  in  the  scan.  There  is  no  doubt  that,  for  "real  time"  appli- 

cations, all  of  those  algorithms  have  a great  disadvantage. 

The  goal  here  is  to  develop  an  algorithm  that  will  detect 
lines,  regardless  of  their  direction. 

2 . The  Algorithm 

The  first,  intuitive,  feeling  is  that  the  algorithm 
of  Section  D,  for  point-detection,  is  correct  for  line 
detection.  But,  if  this  algorithm  is  used,  only  diagonal 
lines  and  the  edge  points  of  vertical  and  horizontal  lines 
can  be  detected.  Therefore,  it  was  necessary  to  improve 
the  algorithm  for  line-detection.  The  problem  lies  in  the 
estimator-equation : 

x(k+l,i+l)  = (k+l,il  + l)  + Gfk+l,ll  + l)  [y  (k+l,i!,  + l)  - (k+1,^  + 1)  ] 

^ correction  term  ^ 

where  the  correction  term  depends  on  the  measurement, 
y(k+l,£+l).  The  "correction  term"  causes  an  error  in  a 
point  where  a target  is  present  (because  in  this  case 
y = X + V + T)  . 

From  another  point  of  view,  looking  on  Fig.  41, 
one  can  see  that  the  basic  assumption  is  that  the  output  of 
the  filter  includes  the  background  alone.  In  that  case,  the 

A 

subtraction  of  the  measured  image,  y(k+l,2.+l),  from  x (the 
output  of  the  filter)  will  result  with  a residual  image 
that  is  T + V.  But  it  was  shown  that  part  of  the  target 
intensity  appears  in  the  output  of  the  filter.  Hence,  after 
subtraction,  the  residual  image  doesn't  include  the  target. 
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This  discussion  indicates  the  improvement  required. 
Recall  that  the  one  step  predictor  is  used  in  the  threshold 
device  to  decide  whether  or  not  a target  is  present.  This 
decision  is  used  to  improve  the  estimation  procedure.  The 
new  algorithm  will  be  a combination  of  the  detection  and 
filtering  process. 


Fig.  46:  The  Combination  Of  Detection  And  Estimation 

For  Improved  Line  Detection 


The  improved  esti.mator  equation  will  be 


x(k+l,i+l)  = 


X (k+1,  ^.+1)  if  r(k+l,Z+l)  > threshold 


x^  ’ (k+l,£+l)  + G(k+1,«.+1)  r(k+l,  +1)  < threshold 

(y(k+l,il+l)  - x'-^' (k+l,il+l)  ] 

5-9 


where,  again 


r(k+l,il+l)  = y(k+l,i+l)  - x (k+1,  Z+1)  . 
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The  idea  behind  this  improvement  is  that  if  the  residual 
is  greate^  than  a certain  threshold  we  know  that  this  is 
due  to  a target,  and  therefore  it  is  wrong  to  include  the 
measurement  y(k+l,£+l)  in  the  estimation  of  that  particular 
point.  In  this  case  the  best  one  can  do  is  to  use  the 
one-step-predictor  as  the  estimator  of  the  point. 

In  the  case  where  the  mean  of  the  target,  T^,  is 
known,  the  equation  for  x(k+l,)l+l)  will  be; 


x(k+l,  2,+D 


f (1) 

x^  ^ (k+l,S.+l) +G(k+1,£+1)  if  r(k+l,2.+l)  > threshold 

• [(y(k+l,£+l)  -T^)  -x^^’  (k+l,Jl+l)] 

1x^^^  (k+l,)l+l) +G(k+l,Jl+l)  if  r(k+l,Jl+l)  < threshold 

• [y(k+l,Jl+l)  -x^^^  (k+l,il+l)] 
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F.  SUMMARY  OF  EXPERIMENTS 
1 . The  Simulated  Image 

In  order  to  check  the  algorithm  of  line  detection, 
a simulated  image  was  created  (the  size;  128  x 128) . 

a)  According  to  the  autocorrelation  function 


R (n,m) 

XX 


^2 


lm| 


the  correlated  background  was  created  by  passing  white  noise 
through  the  filter. 
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x(k+l,«.+l) 


= p^x(k+l,^)  -p^P2X(k,i)  +p2x(k,^+l)  +u(k,£) 


where 

u(k,)l)  is  Gaussian  white  noise,  zero  mean,  and 
variance  of  /Q,  where 

Q = (1  - p^^)  (1  - P2^) . 

b.  A Gaussian  white  noise  with  zero  mean  and 

variance  a (a  = /r)  was  added  to  the  correlated  field, 
n n 

c.  A target  that  has  the  shape: 

N.  P.  G.  S. 

was  added  to  the  image  of  b) . 

d.  The  simulated  image  will  be 

y(k+l,£+l)  = x(k+l,«,+l)  + v(k+l,2,+l)  + T(k+l,^+l) 


Background  white  noise  target 


2 . Detection  Of  The  Target 

a)  The  background  was  estimated  by  the  filter: 

x^^^  (k+1,  H+1)  = p j^x  (k+1,  A)  - Pj^P2X  (k,  ?.)  + P2X  (k,  «,+!) 
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with 


x(k+l,^+l) 


(k+l,Jl+l)  if  r(k+l,^+l)  > TR 

x^^^  (k+l,Jl+l)  + if  r(k+l,^+l)  < TR 

G[y{k+l,il+l)  -x^^^  (k+1,^+1)] 


where : 


r(k+l,«,+l)  = y(k+l,£+l)  - x (k+1 , ^+l)  . 

b)  The  residual  image  r(k+l,Jl+l)  was  passed  through 
a half-wave  rectifier,  and  then  displayed.  The  reason  for 
the  rectification  was  the  fact  that  it  is  known  that 
negative  values  in  the  residual  image  are  only  noise. 

The  display  used  a line-printer  with  letters  representing 
eight  gray  levels. 

c)  The  residual  image  was  fed  into  the  threshold 
device  to  decide  if  there  is  a target  at  point  (k+l,£.+l). 

if  r(k+l,Jl+l)  > TR  — presence  of  target 

if  r(k+l,il+l)  < TR absence  of  target 

This  decision  was  fed  back  to  the  filter  that  estimates  the 
background.  Fig.  47  summarizes  the  procedure  of  simulating 
an  image  and  target  detection. 
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Fig.  47:  Simulation  of  An  Image  and  Detection  of  Targets 


Results  are  given  in  Table  4,  and  Figs.  48-53. 


Table  4 


Case  1 

Case  2 

Case  3 

a 

1 i 

1 

1 

Pi  = P2 

0.96 

0.96 

0.96 

/R  = a 

n 

0.0 

0.2 

0.33 

T 

1.0 

1.0 

1.8 

TR 

0.6 

0.6 

0.9 

G (gain) 

1 

0.46 

0.32 

Results  Are 
Shown  In 
Figures 

43,44 

45,46 

47,48 
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Fig.  48:  Target  Detection  Case  1:  The  Original  Image 
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Fig.  51:  Target  Detection  Case  2:  The  Residual  Image 
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• • • M«* 


v4«  'w  tftio*  d.d«* 


Fig.  53:  Target  Detection  Case  3:  The  Residual  Image 


I 
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VI.  FURTHER  IDEAS  IN  THE  SUBJECT  OF 
RECURSIVE  IMAGE  PROCESSING 

The  purpose  of  this  chapter  is  to  suggest  some 
additional  topics  in  this  research. 

A.  PROBLEMS  IN  RECURSIVE  FILTERING 

It  is  our  belief  that  the  filter  that  was  developed 
in  Chpt.  IV  is  the  best  suboptimal  filter  in  the  structure 
of  (4-12,  4-15).  It  is  so  because  our  approximation  is 
assumed  to  be  correct  in  steady  state.  But  it  is  still  an 
open  question  to  prove  that  this  filter  is  really  the  best. 

In  recursive  estimation  of  more  complicated  structures 
the  "state  space  structure",  as  in  (2-116),  should  be 
developed.  A proposed  filter  might  be: 


(k+l,ji) 

Al 

: ^ 
1 

M(k,l) 

N^^^  (k,£+l) 

• ^ 

N(k,Jl) 

and 


M(k+1,£) 

(k+l,£) 

G^(k,£) 

lU 

(k,£) 

- 

+ 

[y(k,£)  - (Cj^  C2) 

N(k+l,i) 

N^^^  (k,)l+l) 

G^(k,£) 

N^^^  (k,£) 
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or: 


(1) 

(k+l,£+l) 

h : 

M(k,Jl+l) 

^ (1) 

(k+l,«,+l) 

V 

A3  i A^ 

0 

N(k+l,£)j 

and 


M(k+l,i+l) 

^ (1) 

(k+l,Hti) 

G (k+l,Z+l) 
m 

N(k+1,«,+1) 

''(1) 
n'  ' 

(k+l,£+l) 

G (k+l,£+l) 
n 

(k+l,Jl+l) 


(y(k+l,£+l)  - (C^ 


] 


(k+l,^+l) 


^2^ 


6-4 


Following  the  procedure  of  Chpt.  IV  the  gains  G^(k,)l), 

G^(k,2.)  can  be  calculated  using  recursive  equations. 

Recall  that  the  problem  with  the  filter  of  Chpt.  IV  was 

• ' 

V ^ ^ 

the  non-exL stance  of  a "closed  set"  of  recursive  formulas 
for  gain  calculation  (an  approximation  had  to  be  done) . 

The  proposed  filter  that  is  suggested  in  this  chapter  suffers 
from  the  seune  difficulty.  Therefore  it  is  our  belief  that 
formulas  for  the  general  case  cannot  be  found.  Two  dimen- 
sional filters  with  the  proposed  structure  above  should  be 
derived  and  checked  for  specific  cases. 
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The  next  topic  for  additional  research  is  sensitivity 
analysis  of  the  algorithm  to  changes  of  the  field  parameters. 


In  order  to  improve  the  estimation  in  the  case  where  the 
parameters  of  the  field  are  not  known,  one  can  use  an 
adaptive  system  as  follows. 


One  can  apply  the  target  detection  algorithm  on  images 
with  an  autocorrelation  function  from  (2-9)  of  the  form 

Rjjjj(n,m)  = cos  (9^n)  P2  cos  (02^) 

This  can  be  done  only  after  developing  a filter  for  this 
case  (see  proposed  structure  of  the  filter  in  Section  A 
of  this  chapter) . 

The  approach  to  target  detection  in  this  thesis  was  to 
use  a feedback  line  from  the  output  of  the  threshold  device 
back  to  the  filter  (see  Fig.  41) . It  was  found  that  for 
low  values  of  threshold,  this  approach  is  unstable.  The 
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question  of  stability  criteria  for  this  two-dimensional, 
non-linear  operation  remains  unsolved.  Using  a feedback 
from  the  threshold  device  to  the  filter  has  other  disadvan- 
tages. A solution  to  this  problem  might  be  to  estimate 
x(k,il)  without  a feedback  from  the  threshold  device.  The 
residual  image  will  then  be 

r(k+l,Jl+l)  = y(k+l,£+l)  -ax(k+l,Z)  -8x(k,£+l) 
where 


a, 8 have  been  found  by  some  optimization  criteria 

Detection  may  be  improved  by  processing  more  than  one  frame 
or  alternatively  tracking  in  time  domain  from  frame  to 
frame  by  a regular  Kalman  filter. 

An  improvement  may  be  achieved  using  the  algorithm  of 
(5-10),  and  assuming  the  intensity  of  the  target  is  known. 
If  the  intensity  of  the  target  is  not  known,  one  can  use 
this  algorithm  by  adding  a "tracking  loop"  on  the  target 
intensity  in  order  to  estimate  this  intensity  level  on-line 
At  last  all  these  ideas  should  be  checked  on  "real 


life"  images. 


APPENDIX  A 


THE  ONE  STEP  PREDICTOR 

Here  we  compare  between  the  recursive  filter  and  the 
optimal  non-recursive  estimation  (by  direct  use  of  the 
orthogonality  principle) . It  will  be  seen  that  the  one- 
step-predictor  is  not  optimal.  In  order  to  prove  it  a 
specific  case  will  be  shown  and  then  we  will  extend  the 
result. 

Part  1;  A Specific  Case. 

The  recursive  procedure  of  gain  calculation  by  the 
algorithm  of  Chpt.  IV,  leads  to  optimal  estimators  for  the 
first  line  and  column.  That  fact  is  obvious  because  the 
first  line  and  column  could  be  treated  as  a one  dimensional 
case,  using  Kalman  Filtering.  Now  assume  the  measurement 
starts  at  (0,0).  Then,  the  first  place  where  a difficulty 
appears  is  in  the  calculation  of  the  "one-step-predictor" 
of  point  (1,1).  We  shall  calculate  P^^^(l,l)  using  two 
methods:  (1)  optimal  non-recursive  form,  and  (2)  by  using 

the  recursive  form,  in  order  to  compare  the  results. 

Given: 


a » 1 

Pi  ■ p2  “ 0*96 
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y(k,A)  = x(k,«,)  + v(k,£) 


I 


I 

I 

I 

i 


jo  ,n/Oorm^O 

Rw(n,m)  = I 

(r=0.64  ,n=0  and  m = 0 

Find ; The  estimator  of  point  (1,1)  by  using  its  three 
neighbors  (0,0),  (0,1),  (1,0). 


Fig.  43;  The  one  step  predictor.  The  Fig.  shows 
the  pixels  that  take  part  in  the 
one-step-predictor  of  pixel  1,1 


Solution  By  Using  the  Orthogonal  Principle 

By  using  Eq.  C-1  (Appendix  C)  the  next  set  of  equations 

is  obtained. 


XX 


^(0,0)  +R 

R^(1,0) 

R^(0,1) 

^^(-1,0) 

R (0,0) +R 

XX  ' 

^xx<-^'l)  I 

l^(0,l) 

R (1,-1) 

XX  ' 

R^(o,o)  +Ry 

0,0 


0,1 


R^(0,-1) 


(R^(-1,0) 
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Now,  using  Eq.  C-2,  the  error  is  calculated; 


P(l,l) 


optimal 


E{§^(1,1) } 


= E{x^(l,l) 


(5:a  y (p,q)  )x(l,l)  } , 
p » q 


where: 


I 


(p,q)  « {(0,0),  (0,1),  (1,0)} 


1 - - 02p  - 


a^P 


P(l,l) 


optimal 


0.2360 
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Solution  By  Recursive  Filtering 

To  simplify  the  calculations,  the  G's  and  P's 
for  points  (0,0),  (0,1),  (1,0)  will  be  calculated  from  the 
Kalman  filter  equation.  (Recall;  the  results  for  the 
first  line  and  column  by  using  the  algorithm  of  Chpt.  IV 
are  the  same  as  in  the  Kalman  Filtering  results.) 


P(0,0)  = 1 


(Initial  Condition) 


G(0,0)  = 


P^^  (0,0) 

P (0,0)  + R 


0.6098 


P(l,l) 


P (0/  0)R 
P^(0,0)  + R 


= 0.3902 


P^(l,2)  = P^(2,l)  = 0^  P(l,l)  + Q = 0.4380 


G(l,2)  = G(2,l)  » 

P (1,2)  + R 


0.4063 


P(l,2)  » P(2,l)  » 


0.2601 


P (1,2)  + R 


Also  the  calculations  for  the  P-  .(1,1)  and  P.  .(1,1)  can 

U ^ 1 A f 0 

be  done  in  a simplified  way. 

Equivalently  to  Eq.  4-20,  in  the  one-dimensional 
filter,  one  can  write: 


&(k+l)  - F(k+l)p2^(k)  -F(k+l)u(k)  + G (k+1)  v (k+1)  - for  & » 0 
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£(Z+1)  = -F(H+l)u(^)  + G(£+l) v(i+l)  - for  k = 0 

therefore : 

Pq^j^(0,0)  = E{t(0,0) 

= P(0,0)  (l  - G(0,l))p^  = 0.2224 

q(0,0)  = E{6(0,0)'^(1,0)  } 

= P(0,0)  (1 -G(1,0)  )P2  = 0.2224 

Pj^^_l(0,l)  = E{^(0,1)  6’(1,0)} 

= (1 -G(1,0))  (1  -G(0,1))P(0,0)Pj^P2 

» 0.1268 

Now,  the  equation  for  the  one  step  predictor  is: 


P^(l,l)  - p P(0,0)p‘  + Q 


(pj^  -pj^pj  P2) 


/ P(1,0) 

^1,0^°'°) 

Pl^_l(0,l) 

/ Pi ' 

Pi^o(O'O) 

P(0,0) 

-P1P2 

\pi^.l(0,l) 

P(0,1)  / 

1 ^2  / 

+ (1-Pj^^)  (1-P2^) 
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0.2634 


1 f''’ 

P (1,1)  is  not  far  away  from  P d » D optimal  ' anyhow 

different. 

Part  2;  Proof  of  the  General  Case 

First,  prove  that  the  optimal  solution  requires  the 
estimation  error  of  the  predictor  to  be  uncorrelated 
(orthogonal)  to  the  estimators  that  take  part  in  the 
equation  of  the  predictor: 

X (k+l,H+l)  = a^^x(k+l,Jl)  + a2x(k,2,)  + 

and  the  error,  (k+l,2.+l)  is  given  by 

(k+l,£+l)  = x^^^  (k+l,£+l)  - x(k+l,£+l) 

= aj^x(k+l,£)  + a2x(k,£)  + a2x(k,£+l)  - x(k+l,£+l) 

we  want  to  minimize 

E{((S'  ' (k+l,£+l)  ]^]  » E{  (Oj^x(k+l,il)+a2X(k,£)+a2x(k,£+l)-x(k+l,£+l) 
which  requires 

^ E{ (k+l,£+l) ]^}  - 0.  i » 1,2,3. 


198 


We  obtain 


= 0 

1)  , (k,£),  (k,i+l)} 

in  the  one  dimensional  case: 

• x(k)  } = 0 

Now  check  if  this  orthogonality  condition  can  be  satisfied. 
1)  In  the  one  dimensional  case: 

E{£^^Nk+l)  x(k)}  = E{[x^^Nk+l)  -x(k+l)]  x(k)} 

• / 

substituting:  x(k+l)  = px(k)  + u(k) 

E{6^^Nk+l)  x(k)}  = E{[px(k)  - px(k)  - u(k)]  x(k)} 

= E{ [p£(k)  - u (k) ] X (k) } 
assuming  x(k)  is  optimal,  we  know  that: 

E{6-(k)  X (k)  } » 0 
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E{£^^Nk+l,  Jl+1)  x(m,n)} 


(m,n)  = { (k+1, 


\ 


also: 


E{u(k)  x(k)}  = 0 


Therefore: 


(k+1)  x(k)  } = 0 

Conclusion : The  condition  is  satisfied  in  the  case  of  the 

one-dimensional  one-step-predictor. 

2)  In  the  two  dimensional  case: 

Check  if  E{^^^Nk+l,£+l)  x(k,£,+l)}  can  be  made  zero. 

E{e  (k+i/Ji+i)  x(k,jy-i)}  = 

= E{  [ (p^x(k+l,£)-Pj_p2x(k,Jl)+P2x(k,  l+l) ) -x  (k+1,  Jl+1)  ] x (k,  ^+1)  } 

= E{  [ (Pj^x  (k+1,  Z)  -p^P2X  (k,  ^)  +P2X  (k,  ^+1)  ) - (p^x  (k+1,  £) 

- Pj^P2X  (k,i) +P2X  (k,  2.+1) +u  (k,  Jl)  ]x(k,«,+l)  } 

= E{  [Pj_<£,(k+l,i)-Pj_P2  (k,X)+P2^(k,i+l)  ] -x(k,a+l)  } 

Now,  because  x(k+l,Jl),  x(k,)l),  x(k,Jl+l)  are  based  on  different 
sets  of  data,  and  so  are  also  the  corresponding  errors,  the 
only  way  to  make  the  last  expression  zero  for  each  k,  is 
to  make  each  part  of  it  zero. 
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E{§  (k,£+l)  x(k,Jl+l)  } 


0 


E{^(k,^)  x(k,£+l)}  = 0 

E{<S(k+l,Jl)  x(k,)l+l)}  = 0 

But  the  last  expression,  E{<^(k+l,Jl)  x(k-,  £+1)}  cannot 
be  made  zero.  Note  that  ^(k+l,£)  is  a result  that  is 
obtained  from  a different  set  of  data  than  x(k,£+l). 
For  example  one  cannot  require  '^(1,0)  and  x(0,l)  to  be 
orthogonal,  knowing  that: 

x(0,l)  = function  of  y (0 , 0)^  y (0 , 1) 

^(1,0)  = function  of  y (0, 0)^  y (1, 0)  . 

Q.E.D. 
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APPENDIX  B 


Let  us  show  how  one  can  write  nonrecursive  equations  for 
X (k+1 , )l+l)  , y(k+l,Jl+l),  x(k+l,H+l)  and  (k+l,J,+l)  and  then 
show  some  properties  of  the  correlation  between  these  values. 
From  the  recursive  equation: 


x(k+l,Jl  + l)  = p^x(k+l,Jl)  - pj^P2x(k,Jl)  +P2x(k,2.+1)  +u(k,i) 


one  can  write: 


initial  conditions 


x(k+l,£+l)  = Za(0,£)x(0,^)  + i:a(k,0)x(k,0] 


+ I a(i,j)u(i,j) 

(i,  j)eJ^2 


= M + 


(i,  j)€i^2 


a(i,  j)u(i,  j) 


B-1 


Now,  using  Eq.  2: 


y(k+l,)l+l)  = 


(i,j)cn2(k,Jl) 


a(i,  j)u(i,  j)  ] + 


v(k+l,^+l) 


B-2 


Using  2,  12,  15,  A-1,  A-2; 
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4 


x(k+l,Jl+l)  = N + E S(i,j)u(i,j) 

(i, j ) £^2 


+ E Y(i/j)v(i,j)  B-3 

(i,  j ) £^2  (k+1, 

Theorem:  E{u  (k,  2.)€.'^  (k,  2.)  } = E{^(k,  X, ) u (k,  X, ) } = 0 

Proof t Using  B~l/  B-3  we  see  that  §,(k,l)  does  not 

include  the  value  u(k,il).  Therefore  by  using  V-ywe  have 
completed  the  proof. 

Theorem:  E{v  (k+1,  (k,  X-)  } = E{& (k,  £ ) v (k+1 , 2,+D  } = 0 

I 

I 

Proof : The  proof  is  in  the  same  procedure  as  above. 

By  substituting  B-1,  B-3  into  the  expression  for  §{k,J,) 
we  see  that  <£lk,H)  does  not  include  the  value  v(k+l,il  + i), 
and  therefore,  by  using  Eq.y-7,  we  have  completed  the  proof. 

In  the  same  way  one  can  show  that: 

E{f(k,£)  u(k+i,)l+j)}  = 0 

E{fe(k,2.)  V (k+i+1,  H+j+1) 

if  i > 0 

... 
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APPENDIX  C 


NON-RECURSIVE  ESTIMATION 

1.  Motivation  For  The  Non-Recursive  Estimation 

The  goal  of  this  thesis  is  recursive  estimation  in 
image  processing.  Chpt.  IV  derived  a recursive  "least 
square  error"  estimator  that  was  not  optimal  (the  error  was 
not  orthogonal  to  the  data) . So,  the  motivation  to 
this  section  is  to  compare  the  recursive  filter  to  the 
optimal  solution. 

2 . The  Basic  Procedure 

The  non  recursive  estimation  is  based  on  the  orthogonality 
principle,  which  is  given  in  Eq.  3-10: 


E{[x(k,i)  - Z a y(p,q)]  y(i,j)} 

(p,q)enj^(k,il) 


0 


C-1 


and  writing  it  in  the  form: 


E{[x(k,^)  - E a (x (o,q)  + V (p,q) ) ] 

(p,q)efi^(k,£) 

• [x(p,q)  + v(p,q) ] } = 0 


V is  the  noise  that  is  added  to  the  state  x. 


I 
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Let  us  clarify  the  procedure  with  an  example; 


Given:  1) 


R (n,m) 

XX  ' ' 


= p. 


n| 


m 


2)  R^(n,m) 


j^O  ifn^O  or  mj^O 

R if  n = 0 and  m = 0 


v(k,il)  is  uncorrelated  to  any  of  the  random- 
field  values,  x(i,j). 

3)  y(k,Z)  = x(k,Z)  + v{k,^) 


Find:  the  four  optimal  coefficients  to  estimate  the  value 

x(0,0)  by  the  set  of  measurements  in  a 2x2  area: 
y(0,0),  y(0,l),  y(l,0),  y(l,l) 


)/(0,o) 

y/"o,i) 

/0,o, 

/r>.0 

yco,.)  ^ 0)  ^ 

Fig.  55:  Non  Recursive  Estimation 

Solution:  If  the  noise  is  uncorrelated  to  the  signal 

x(k,4),  then  the  set  of  equations  that  follows 
from  the  orthogonality  condition  is: 
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4l(-l,-l)\  /r^(0,0)+R  R^(1,0)  ' 

/^xx^"^'°^  R^(0,0)+R  ^xx^°'"^^ 

R(-1,0)  Rj^(0,0)+R  Rjoj(l,0) 

,R(0,0)/  \ Rxx^"^''^>  ^xx^°'"^^  ^xx^"^'°>  Rjoj(0'0)+R^ 


/ 1+R 


P^P2 


P1P2 


P^P2 


P1P2 


From  this  four  equation  system,  one  can  find  the  a's. 
3.  The  Estimation  Error  (variance) 


z(k,l)  = x(k,^)  - X = x(k,Z)  - Z a y(p,q) 


e {k,l)  = [x(k,£.)  - Z(x  yCp^q)] 

P r4 


E{e^}  * E{x^(k,4)  - 2x(k,*,)i:a  _y(p,q) 

P » 4 


+ {la  v(p,q))^} 

P f 4 


206 


E{  = E{x^(k,)l)  -x(k,?,)Ea  y(p»q) 

p » q 


+ a_  „y(p»q)(J:a  _y(p,q)  - x)} 

H ^ 4 P f 4 


The  zero  was  set  due  to  the  orthogonality  principle. 


P(k,l)  = E{  ^(k,^)}  = E{x^(k,i)  - (Za  y (p,q) ) x (k, i ) } 

p / q 


Exeunple ; If  we  continue  the  previous  example  (the  estimation 
of  a point  by  a 2x2  area)  : 


E{  } - 1 - «o,0^1^2  ■ “0, 1^2  " “l, 0*^1  ‘°i,l 


Theorem;  If  the  measurement  y(k,J,)  is  included  in  the 
estimation  of  the  point  x(k,)l),  then: 


P(k,i)  = E{  ^(k,i)}  = R 


Proof ; Assume  the  estimation  of  the  point  x(k,2.)  is  done 
without  the  point  Ck,l},  by  using  a set  of  measurements  Y. 

A 

Let's  call  this  estimator  x (p  = Prediction). 

P 


Xp(k,z) 


^ a.  ^y(i/j) 

(i,j)eY 


(i,j)  {k,^) 
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Xp(k,£)  is  the  "predictor"  of  x(k,il). 

Also  assume  that  the  variance  of  the  error  of  the  predictor 

is  P^.  Now  assume  another  estimator  that  Uses  the 

set  y and  the  measurement  y^k,^).  in  that  case  we  can  say; 

x(k,A)  = Xp(k,X,)  + y(k,£). 


x(k,2.)  is  supposed  to  be  an  optimal  combination  of 
x(k,£)  and  y(k,2,)  . Because  there  is  no  correlation 

between  v(k,£)  and  other  values,  this  is  a case  of  a com- 
bination between  two  estimators  that  have  uncorrelated 
errors  [see  Ref.  4,  Chapter  on  optimal  smoothing]. 


estimator  A: 


Xp(k,«,) 


with  variance  P . 

P 


estimator  B;  y(k,2.)  with  variance  R. 


The  optimal  combination  is  [due  to  Ref.  4] : 


x(k,)l) 


R^rp-  + rtV 

P P 


with  variance; 


P(k,^) 


P R 
„ P ^ 

ft  + P 


P 


Q • £)  • D • 
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4.  Unbiased  Estimator 


Given ; A random  field  x(k,2.)  with  mean  x. 

A noise  measurement  y{k,l)  = x(k,«.)  + v(k,)l). 

We  wish  to  find  an  unbaised  estimator  for  this  case. 

Define : x(k,J.)  - x'(k,Jl)  + x 

x*(k,J,)  is  a random  field  with  mean  zero.  Suppose 
we  want  to  estimate  x'(k,)l).  To  do  that  define 

y ’(k,Jl)  ^ x'(k,i)  + v(k,i) 

= y(k,Z)  - X 


the  estimation  of  x'(k,2.)  is; 

X (k, H)  = lay* (p,q) 

P ^ 4 

x(k,£)  is  unbiased  because  x'(k,)l)  is  unbiased. 
And  now,  the  estimation  of  x(k,)l)  is: 

x(k,£)  = x’(k,£)  + X = y'(Pfq)  + X 

P f 4 


a - 

P / 4 


x(k,£)  = Ea  _ y_  _ 

P»q  Pf<3 


x)  ] + X 

+ x[l  - Ea„  ] 
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Conclusion;  The  estimator  can  be  done  unbaised  regardless 
what  the  coefficients  a are.  Therefore,  the  calcxalation 

p»q 

of  the  a's  is  the  same  as  before.  One  has  only  to  add  one 
term  to  the  estimator  equation.  This  term  depends  on  the 
mean. 
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